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Abstract 


A  surface  panel  method  suitable  for  the  analysis  of  marine  propellers  is  developed 
and  applied  to  various  geometries  to  demonstrate  its  e:Tectiveness.  A  reliable  pressure 
distribution  around  a  marine  propeller,  especially  near  the  leading  edge  of  the  blade, 
is  obtained.  Hub  effects  are  naturally  included  by  distributing  panels  on  the  hub 
surface. 

Detailed  study  of  the  flows  at  the  trailing  edge  near  the  tip  suggests  a  pressure 
Kutta  condition,  which  requires  the  pressures  of  the  last  panels  at  the  trailing  edge 
be  equal.  Due  to  the  nonlinear  aspect  of  the  pressure  Kutta  condition,  an  iterative 
process  is  employed.  An  efficient  approximation  of  the  ultimate  wake  is  achieved  by 
replacing  it  with  a  sink  disk  at  the  beginning  of  the  ultimate  wake. 

The  sample  geometries  include  an  ellipsoid  at  zero  angle  of  attack,  a  circular 
planform  wing,  a  rectangular  planform  wing  with  varying  sw'eep  angles,  a  wing-body 
configuration,  a  long  axisymmetric  duct,  and  a  marine  propeller.  Calculated  pressure 
distributions  around  the  wing-body  configuration  are  in  excellent  agreement  with  the 
experimental  data.  Calculated  thrust  and  torque  for  the  propeller  agree  well  with 
experimental  results. 
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Chapter  1 
Introduction 


Calculation  of  the  pressure  distribution  around  a  marine  propeller  is  a  fundamental 
goal  of  many  naval  architects,  yet  a  successful  calculation  of  a  reliable  pressure  dis¬ 
tribution,  especially  near  the  leading  edge  and  the  blade  tip,  has  not  been  reported. 
Not  only  should  the  thrust  and  the  torque,  which  are  calculated  by  the  integration 
of  the  pressure  distribution  on  the  blade  surface,  be  accurate,  but  the  local  minimum 
pressure  at  the  leading  edge  of  the  blade  should  also  be  accurate  to  calculate  cavita¬ 
tion  inception  and  extent.  \Ioreover,  the  locations  of  the  local  stagnation  point  and 
the  minimum  pressure  point  at  the  leading  edge  of  the  blade  are  essential  inputs  to 
the  calculation  of  the  boundary  layer  on  the  blade  surface. 

Most  of  the  design  and  analysis  of  marine  propellers  rely  heavily  on  lifting  surface 
theory.  Many  lifting  surface  codes,  which  use  a  discrete  vortex/source  line  represen¬ 
tation  of  the  blade  on  the  mean  camber  surface,  have  been  developed  and  used  in 
the  prediction  of  the  steady/unsteady  performance  of  marine  propellers  [14), [5].  As 
a  consequence  of  the  linear  superposition  of  the  thickness  and  the  lifting  problem  in 
lifting  surface  theory,  the  prediction  of  the  pressure  distribution  at  the  leading  edge 
is  not  valid.  Lighthill’s  correction  may  be  applied  to  the  local  flow  around  the  leading 
edge,  however  the  validity  is  questionable  for  general  three-dimensional  flows.  More¬ 
over  the  complete  exclusion  of  the  hub  in  the  numerical  model  of  the  lifting  surface 
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theory  has  been  questioned.  Recent  research  on  the  effect  of  the  hub  on  the  perfor¬ 
mance  of  a  propeller  by  Wang  [23]  showed  that  the  effect  could  be  significant  in  those 
propellers  with  a  large  hub  radius  or  a  small  number  of  blades.  These  difficulties  can 
be  resolved  naturally  by  adopting  a  surface  panel  method  and  distributing  panels  on 
the  exact  propeller  surface  including  the  hub  surface.  This  becomes  possible  as  high 
speed  computers,  with  large  capacity,  become  more  commonly  available. 

Many  kinds  of  panel  methods  have  been  in  use  for  aerodynamic/hydrodynamic 
applications  since  Hess  [8]  proposed  the  surface  source  method.  Although  all  the 
properly  formulated  panel  methods  are  exact  in  the  sense  that  the  numerical  solutions 
converge  to  the  common  solution  as  the  number  of  panels  is  increased,  this  does  not 
imply  that  all  the  panel  methods  are  equally  successful.  Indeed,  vast  differences  exist 
with  respect  to  the  prediction  accuracy  versus  computational  effort  and  the  reliability 
for  extreme  geometries. 

In  the  present  paper,  the  characteristics  of  the  various  panel  methods  are  reviewed 
and  compared.  As  a  result,  a  low  order  panel  method  bcised  on  the  perturbation 
potential  is  chosen  because  of  its  robustness  with  respect  to  extreme  geometries,  and 
its  relatively  smaller  computational  effort. 

Extreme  geometry  of  the  blade  of  marine  propeller  had  prevented  the  surface  panel 
method  from  being  applied  to  the  marine  propellers.  To  the  author’s  knowledge,  only 
two  attempts  have  been  made  to  apply  the  panel  method.  Hess  [9]  extended  his 
surface  source  method  to  this  problem.  But  the  surface  source  method  suffers  as  the 
thickness  of  the  blade  becomes  small,  and  the  results  toward  the  thin  blade  tip  become 
questionable.  Koyama  et.  al.  [15]  applied  a  Morino  type  low  order  panel  method 
to  marine  propellers.  But  his  pressure  distribution  near  the  tip  showed  a  spurious 
positive  loading  at  the  trailing  edge.  This  is  due  to  the  failure  of  the  Morino ’s  Kutta 
condition  which  does  not  account  for  the  three  dimensional  cross  flow  effects.  In  the 
present  work,  extensive  attention  is  given  to  the  three  dimensional  trailing  edge  flow. 
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As  a  result,  a  pressure  Kutta  condition,  which  requires  the  pressures  of  the  last  panels 
at  the  trailing  edge  be  equal,  is  suggested.  Due  to  the  nonlinear  aspect  of  the  pressure 
Kutta  condition,  an  iterative  solution  procedure  is  employed. 

Numerical  calculation  is  performed  for  various  geometries  using  the  selected  panel 
method  with  the  pressure  Kutta  condition.  The  examples  selected  in  the  present 
work  include  an  ellipsoid  at  zero  angle  of  attack,  a  circular  planform  wing  with  vary¬ 
ing  thickness,  a  rectangular  planform  wing  with  different  sweep  angles,  a  wing-body 
configuration,  a  long  axisymmetric  duct,  and  a  marine  propeller. 
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Chapter  2 


Fundamentals  of  the  panel  methods 

2.1  Introduction 

Panel  methods  have  been  in  use  for  certain  aerodynamic/hydrodynamic  applications 
since  Hess  [8]  proposed  the  surface  source  method  in  1964.  Since  then,  various  other 
formulations  have  appeared  which  offer  advantages  in  terms  of  accuracy,  computa¬ 
tional  efficiency  or  versatility.  Although  all  the  properly  formulated  panel  methods 
are  exact  in  the  sense  that  the  numerical  solutions  converge  to  the  common  solution 
as  the  number  of  panels  is  increased,  this  does  not  imply  that  all  the  panel  methods 
are  equally  successful.  Indeed,  vast  differences  exist  with  respect  to  the  prediction 
accuracy  versus  computational  effort,  reliability  and  simplicity. 

In  this  chapter,  the  basic  mathematical  theory  behind  the  various  panel  methods 
is  reviewed  and  the  characteristics  of  each  method  are  compared  in  order  to  determine 
the  most  suitable  method  for  the  analysis  of  marine  propellers.  The  common  basis  of 
the  apparently  different  panel  methods  will  then  be  evident,  suggesting  the  following 
grouping  of  the  panel  methods. 

•  Potential  field  formulation. 

-  Perturbation  potential  method. 

—  Total  potential  method. 
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•  Velocity  field  formulation. 


-  Mixed  source  and  dipole  method 

-  Dipole  method.  (Equivalently,  vorticity  method) 

-  Source  based  method. 


2.2  Statement  of  the  problem 


Consider  a  closed  three  dimensional  domain  V  with  boundary  S,  the  unit  normal 


vector  n  to  5  being  oriented  into  V,  as  shown  in  Figure  2.1.  The  boundary  S  is 
composed  of  the  body  surface  Sb,  the  wake  surface  Sw,  and  the  outer  control  sur¬ 
face  S(x)  surrounding  the  body  and  wake  surface.  The  body  is  subject  to  the  inflow 
velocity  Uoo.  With  the  assumptions  that  the  fluid  in  V  is  incompressible,  inviscid, 
and  irrotational,  there  exists  a  perturbation  velocity  potential  4>  which  satisfies  the 
Laplace  equation. 


vV  =  0. 


(2.1) 


A  boundary  value  problem  can  be  constructed  by  specifying  boundary  conditions 
on  the  boundary  S  as  follows: 

•  The  kinematic  boundary  condition  should  be  satisfied  on  the  solid  body  surface 


Sb, 


(2.2) 


•  The  wake  surface  Sw  is  assumed  to  have  zero  thickness.  The  normal  velocity 
jump  and  the  pressure  jump  across  Sw  is  zero,  while  a  jump  in  the  potential  is 
allowed. 


(^P)on  Sw  =  -P  =  0, 


(2.3) 

(2.4) 
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Figure  2.1:  Notation  for  a  general  body  for  the  application  of  Green’s  theorem. 
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For  the  steady  lifting  problem,  the  potential  jump  across  the  wake  surface  is 
the  same  as  the  circulation  around  the  body,  and  is  constant  in  the  streamwise 
direction  on 

{A4>)onS,=<i>^  -r  =T.  (2.5) 

•  A  Kutta  condition  is  required  at  the  trailing  edge  to  uniquely  specify  the  circu¬ 
lation.  In  its  most  general  form,  it  states  that  the  flow  velocity  at  the  trailing 
edge  remains  bounded:  i.e., 

<  oo.  (2.6) 

•  On  the  outer  control  surface  S^o,  the  perturbation  velocity  due  to  the  body 
should  vanish  in  the  limit  where  this  surface  is  an  infinite  distance  from  the 
body. 

V(^  ->  0,  as  Soo  oo  (2.7) 


According  to  Lamb  [16],  this  boundary  value  problem  for  the  velocity  potential 
outside  the  body  surface  can  be  transformed  into  an  integral  equation,  upon  consid¬ 
eration  of  a  fictitious  fluid  in  V,  which  is  the  domain  internal  to  the  body  surface  Sb- 
Thus  for  the  field  point  p  in  V 


47r<^(p)  =  11  -  4>'{<l)) 


where 


Sb 

+ 


1 1  A4>{q) 

Svv 


dn,  i?(p;  q) 

1 


-( 


d<i>{q)  _  d(p'{q)  _ ]_ 

an,  an,  i?(p;g) 


dS 


dn,R{p;  q) 


dS, 


(2.8) 


d>  =  perturbation  velocity  potential  in  V , 

4>'  =  perturbation  velocity  potential  in  V', 

p(x,y,2)  =  field  point  where  induced  potential  is  calculated, 
9(^)^)f)  =  source  point  where  singularity  is  located, 

R[p\q)  —  distance  between  point  p  and  q, 

=  ^/(x  -  0^  +  (y  -  pY  +  {z-  c)^ 
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^  =  normal  derivative  with  respect  to  the  point  q. 

This  equation  may  be  regarded  as  a  representation  of  the  velocity  potential  in 
terms  of  a  normal  dipole  distribution  of  strength  [(f>  —  <j)')  on  the  body  surface  Sb,  a 
source  distribution  of  strength  (f^  -  |^)  on  5b,  and  a  normal  dipole  distribution  of 
strength  A<^  on  the  wake  surface  S\y. 

Since  the  fictitious  fluid  inside  Sb  does  not  have  physical  meaning,  we  can  choose 
the  internal  velocity  potential  4>'  to  suit  our  convenience.  Thus,  by  choosing  an  appro¬ 
priate  <f>^  in  Equation  2.8,  we  can  formulate  various  panel  methods  which  use  different 
sets  of  singularities, 

2.3  Potential  field  formulation 


If  we  choose  the  fictitious  potential  as  0'  =  0  on  Sb,  Equation  2.8  for  the  field  point 
p  on  the  body  surface  Sb  becomes 


2n4>{p) 


11  <^(9) 

Sb 


d  1 
5n,  R{p\q) 


d(b{q)  1 
an,  R[p-,q) 


+Jj 

Sw 


d  1 

an,  i?(p;  q) 


dS. 


(2.9) 


Here  the  surface  integral  on  Sb  must  be  defined  to  exclude  the  immediate  vicinity  of 
the  singular  point.  By  choosing  the  internal  potential  on  Sb  cj)'  =  0,  the  internal 
flow  can  be  shown  to  be  an  undisturbed  flow  of  =  (f>^,  where  the  total  internal 
velocity  potential  is  defined  as  $  '  =  d>oo  +  <t>'. 

Because  is  known  on  Sb  from  the  boundary  condition  (Equation  2.2),  Equa¬ 
tion  2.9  is  a  Fredholm  integral  equation  of  the  second  kind  for  the  dipole  strength  4), 
which  is  also  the  potential  value  on  the  body  surface  Sb-  The  potential  jump  across 
the  wake  surface  can  be  set  equal  to  the  difference  between  the  potential  values  of 
the  upper  and  lower  surfaces  at  the  trailing  edge,  which  replaces  the  Kutta  condi¬ 
tion.  Discretization  of  Equation  2.9  will  lead  to  a  linear  system  of  equations  for  the 
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unknown  <f>.  The  surface  velocity,  hence  the  pressure,  on  Sb  can  be  calculated  by  a 
numerical  differentiation  of  the  potential  distribution.  This  form  of  panel  method  was 
introduced  by  Morino  [20],  and  adopted  in  the  present  paper.  We  will  refer  to  this  as 
Morino’s  method  or  the  perturbation  potential  method. 

If  there  exists  an  inflow  velocity  potential,  4>oo,  such  that  Vd>oo  =  we  can 
formulate  another  form  of  the  panel  method  by  choosing  the  internal  potential  in 
Equation  2.8  cis  a  negative  of  the  inflow  velocity  potential,  i.e.,  <f)'  =  The 

source  strength  in  Equation  2.8  becomes  zero  because  of  the  boundary  condition 
(Equation  2.2),  while  the  dipole  strength,  which  is  the  difference  between  <j)  and  4>', 
becomes  the  total  potential, 

4>  —  <p'  =  (f>  +  (j>oo  =  (2-10) 

As  the  point  p  approaches  the  body  surface  Sb,  the  contribution  from  the  imme¬ 
diate  surface  on  Sb  in  the  first  term  of  Equation  2.8  is 

J;”„  .li-?.  //(«'  - 

The  resulting  equation  is 

Sb 

Sw 

This  equation  can  be  regarded  as  a  representation  of  the  total  velocity  potential 
in  terms  of  a  normal  dipole  distribution  only  on  the  body  surface  Sb  and  the  wake 
surface  Sw-  Given  the  inflow  velocity  potential  values,  this  is  also  a  Fredholm  integral 
equation  of  the  second  kind  for  the  total  potential  $.  Discretization  of  this  equation 
gives  another  form  of  panel  method,  which  we  will  refer  to  ais  the  total  potential 
method. 
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2.4  Velocity  Field  Formulation 


Instead  of  forming  an  integral  equation  in  the  potential  field,  we  can  alternatively 
construct  one  in  the  velocity  field.  Taking  the  normal  derivative  of  Equation  2.8  with 
respect  to  the  field  point  p,  the  resulting  equation,  when  the  field  point  p  is  on  Sb,  is 


47r 


dn„ 


-  2^ a(p)  +  j j  a(„)-^(  j^)dS  +  j I  ^U) 

Sb  ^  ’  Sb 


dnpdn^  R[p\  q) 


+ 


Jl  M{q) 

Sw 


dupdn,  R[p\  q) 


dS, 


dS 


(2.13) 


where  ^  ^  ~  ^  and  p  =  <f>  -  4> . 

It  can  be  shown  that  a  dipole  distribution  on  a  closed  body  surface  is  equivalent  to 
a  vorticity  distribution  with  strength  7,  which  is  calculated  as  a  vector  product  of  the 
local  surface  gradient  of  the  dipole  strength  and  the  normal  vector,  (see  Appendix 
B)  Thus  we  can  alternatively  write  Equation  2,13  as 

=  S’' +  //''(?)  +  9) 

Sb  Sb 

+  11  np-  7(g)  X  Vp-^d5,  (2.14) 

Sw 


where  7  =  n,  x  V2_<j(d>  -  <i)')  . 

Here  again,  by  choosing  different  values  for  the  internal  potential  4>' ,  we  can  express 
the  normcd  velocity  on  Sb  in  terms  of  different  sets  of  singularities. 

If  we  choose  the  internal  potential  in  Equation  2.13  as  (;!»'  =  0,  it  can  be  shown 
that  =0  on  Sb.  Then, 


27r 


d<;!>(p) 

dn„ 


pdn^  R 


+  If  — ;5 — 

J  J  OTlpOTlg  R 

Sw 


(2.15) 
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where  o’  =  H  =■  4>. 

The  right  hand  side  of  the  Equation  2.15  can  be  regarded  cis  the  normal  induced 
velocity  at  p  due  to  the  mixed  distribution  of  normal  dipoles  of  strength  4>  on  Sb, 
sources  of  strength  on  5s,  and  normal  dipoles  of  strength  A4>  on  5^.  The  source 
strength  and  the  left  hand  side  of  Equation  2.15  is  given  by  the  kinematic  boundary 
condition  on  the  body  surface.  Thus  Equation  2.15  is  an  integral  equation  of  the  first 
kind  for  the  unknown  dipole  strength  p. 

Discretization  of  Equation  2.15  gives  another  form  of  panel  method,  which  we  will 
refer  to  zs  the  mixed  source  and  dipole  method.  This  can  be  regarded  cls  the  velocity 
field  formulation  of  the  perturbation  potential  method. 

An  equivalent  formulation  derived  from  Equation  2.14  leads  to  a  mixed  distribution 
of  sources  and  vortices  instead  of  dipoles. 


27r 


d4>{p) 

dup 


Sb  ^  Sb 

+  // • '7(9)  X  Vp-^<i5,  (2.16) 

Sw 


where  cr  =  ^ind  ^  =  rig  x 

If  we  choose  the  internal  potential  on  5s  as  then  <^'  can  be  shown  to  be 

equal  to  — 0oo  throughout  the  inside  of  the  body.  Then  the  source  strength  becomes 
zero  and  Equation  2.13  becomes 


Sb  Sw 


where  /x  =  (^  +  <^00  = 

This  is  also  an  integral  equation  of  the  first  kind  for  the  unknown  dipole  strength 
p.  The  panel  method  derived  from  Equation  2.17  will  be  referred  to  as  the  dipole 
method.  This  can  be  regarded  as  the  velocity  field  formulation  of  the  total  potential 
method. 


22 


Due  to  the  equivalence  between  dipoles  and  vortices,  Equation  2.17  can  be  written 
in  a  different  form, 

4'^=  //  (2.18) 

Sb  5h. 

where  7  =  n,  x  V2_^$.  We  will  refer  to  the  panel  method  derived  from  this  equation 
cLS  the  vorticity  method. 

If  we  choose  the  vorticity  strength  in  Equation  2.14,  such  that  it  has  a  given  shape 

function  g(t)  along  the  chordwise  panels  and  the  spanwise  circulation  r(s)  is  yet  to 

be  determined,  then 

"  Sb  ^  Sb 

+  // -7  X  Vp— d5,  (2.19) 

Sw 

where  7  =  r(s)3(t)t.,,  t-,  is  the  direction  of  the  vorticity,  and  s  and  t  are  the  spanwise 
and  chordwise  coordinates.  Given  the  vorticity  shape  function  g{t)  and  the  normal 
velocity  on  Sb  from  the  boundary  condition,  this  is  a  Fredholm  integral  equation 
of  the  second  kind  for  the  unknown  source  strength  a  and  the  spanwise  circulation 
distribution  r(s).  This  is  the  form  of  the  original  surface  source  method  by  Hess  '8]. 
We  will  refer  to  this  form  as  the  source  based  method. 

2.5  Comparison  of  the  characteristics  of  the  vari¬ 
ous  panel  methods 

To  compare  the  characteristics  of  each  panel  method,  the  above  described  methods  are 
implemented  as  numerical  codes  for  the  analysis  of  the  two-dimensional  flow  around 
a  hydrofoil.  It  is  assumed  that  the  characteristics  of  each  panel  method  are  preserved 
for  the  two-dimensional  case,  even  though  the  principal  objective  of  the  present  study 
is  to  compare  those  for  the  three  dimensional  case. 
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Numerical  implementation  of  the  panel  methods  involves  approximations  of  vari¬ 
ous  sorts  and  inevitably  introduces  discretization  errors  which  would  decreaise  as  the 
number  of  panels  is  increased.  The  principal  approximations  involved  with  any  panel 
methods  can  be  summarized  as  follows: 

•  Discretization  of  the  geometry, 

•  Discretization  of  the  singularity  distributions, 

•  Boundary  conditions  are  satisfied  at  discrete  collocation  points, 

•  A  numerical  Kutta  condition  must  be  imposed. 

The  simplest  discretization  will  be  used  in  the  numerical  implementation  of  the  panel 
methods  in  order  to  compare  the  characteristics  of  each  method. 

Of  the  five  different  panel  methods  which  were  described  in  the  previous  sections, 
the  source  based  method  has  been  used  most  widely  and  is  known  to  be  accurate  for 
most  geometries.  However,  the  source  based  method  is  also  known  to  be  inaccurate 
for  relatively  thin  foil  sections.  For  marine  propeller  applications  this  is  a  severe 
disadvantage  because  a  typical  thickness/chord  ratio  is  as  little  as  two  percent  near 
the  tip.  Moreover,  the  source  based  method  is  known  to  be  inaccurate  for  internal 
flows,  such  as  the  flow  inside  a  long  duct. 

Because  the  characteristics  of  the  source  based  panel  method  are  relatively  well 
known,  four  other  panel  methods  are  actually  implemented.  They  are  the  perturbation 
potential  method,  the  total  potential  method,  the  mixed  vortex  and  source  method, 
and  the  vortex  method.  A  brief  description  of  each  method  will  be  given,  and  the 
more  detailed  numerical  implementation  is  included  in  Appendix  A. 

In  the  perturbation  potential  method,  the  foil  geometry  is  replaced  by  an  N-faced 
polygon,  where  N  is  the  number  of  panels.  The  singularity  strength  on  each  panel 
is  assumed  to  be  piecewise  constant.  The  collocation  point,  where  the  discretized 
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integral  equation  is  satisfied,  is  selected  as  a  midpoint  of  each  panel.  A  numerical 
Kutta  condition  is  satisfied  by  requiring  that  the  potential  jump  in  the  wake  be  the 
difference  between  the  potential  values  of  the  upper  and  lower  trailing  edge  panels. 
The  above  treatment  results  in  a  system  of  linear  equations  for  the  unknown  dipole 
strengths. 

The  surface  velocity,  hence  the  pressure,  is  calculated  by  a  second  order  differen¬ 
tiation  of  the  resulting  perturbation  potential.  Finally,  the  lift  and  drag  are  obtained 
by  summing  the  element  pressure  forces  which  are  calculated  by  multiplying  the  pres¬ 
sure  at  the  midpoint  by  the  arc-length  of  each  panel.  Alternatively,  the  lift  can  be 
calculated  from  Kutta-Joukowsky’s  law  as  pUT ,  where  F  is  given  as  the  potential 
jump  on  the  wake  surface,  and  the  drag  should  be  zero. 

The  computer  code  is  applied  for  a  typical  hydrofoil  section  with  a  thickness/chord 
ratio  of  four  percent  and  a  camber/chord  ratio  of  two  percent.  The  thickness  form 
is  chosen  as  the  NACA  66  mod.  form  and  the  camber  distribution  is  chosen  as  the 
a=0.8  mean  camber  line,  which  is  widely  used  in  marine  propeller  blades  because  of 
its  good  performance  with  respect  to  cavitation  inception.  The  panel  arrangement 
for  40  elements  is  illustated  in  Figure  2.2. 

Convergence  characteristics  of  the  perturbation  velocity  potential  distribution  and 
the  pressure  distribution  obtained  by  the  perturbation  potential  method  are  shown  in 
Figure  2.3  and  2.4,  where  number  of  panels  is  increased  as  20,  40,  80  and  160.  The 
results  with  40  panels,  which  is  a  typical  number  of  the  chordwise  panels  in  three- 
dimensional  application,  are  shown  to  be  very  close  to  those  with  160  panels.  The 
pressure  distribution  with  160  panels  is  considered  converged  and  will  be  used  as  a 
standard  datum  for  comparison  with  the  other  panel  methods. 

Geometric  discretization  for  the  total  potential  method  is  the  same  as  that  for  the 
perturbation  potential  method.  Instead  of  distributing  sources  and  normal  dipoles, 
only  normal  dipoles  are  distributed  on  the  foil  surface.  Since  the  kernel  of  Equa- 
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Figure  2.2:  Panel  arrangement  of  a  typical  foil  section.  (NACA  66  mod.  +  a=0.8 
mean  camber  line,  t/c=0.04,  f/c=0.02) 
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Figure  2.3:  Convergence  characteristics  of  the  perturbation  potential  distribution 
calculated  by  the  perturbation  potential  method.  (NACA  66  +  a=0.8,  t/c=0.04, 
f/c=0.02,  q=1.5  deg) 
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Figure  2.4:  Convergence  characteristics  of  the  pressure  distribution  calculated  by  the 
perturbation  potential  method.  (NACA  66  -f  a=0.8,  t/c=0.04,  f/c=0.02,  a=1.5  deg) 
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tion  2.12  is  the  same  as  that  of  the  Equation  2.9,  the  influence  coefficient  matrix, 
which  is  the  left-hand  side  of  the  linear  system  of  equations,  is  identical  to  that  of 
the  perturbation  potential  method.  Only  the  right-hand  side  of  the  linear  system  is 
different  and  the  unknown  is  the  total  potential  rather  than  the  perturbation  poten¬ 
tial.  The  potential  jump  on  the  wake  surface  is  equal  to  the  difference  of  the  total 
potential  values  of  the  upper  and  lower  trailing  edge  panels. 

Convergence  characteristics  of  the  total  velocity  potential  distribution  and  the 
pressure  distribution  by  the  total  potential  method  are  shown  in  Figures  2.5  and  2.6. 
Comparing  Figures  2.3  and  2.5  we  notice  that  the  convergence  of  the  total  potential 
method  is  much  faster  than  that  of  the  perturbation  potential  results.  However,  the 
potential  jump  in  the  wake,  or  the  circulation  around  the  foil,  by  the  total  potential 
method  converges  at  the  same  rate  as  that  of  the  perturbation  potential  method. 
The  fast  convergence  of  the  total  potential  in  Figure  2.5  is  attributable  to  the  fact 
that  the  inflow  potential  is  included  in  the  total  potential.  The  disadvantage  of  the 
total  potential  method  is  that  for  a  non-uniform  inflow  velocity,  the  inflow  velocity 
potential  is  either  not  defined  or  difficult  to  define.  In  such  a  case  this  method  can 
not  be  applied. 

Table  2.1  provides  a  comparison  of  the  computed  values  of  lift  and  drag  coeffi¬ 
cients  obtained  by  an  integration  of  the  pressure  distribution,  with  the  lift  coefficient 
calculated  by  Kutta-Joukowsky’s  law. 

The  velocity  field  panel  methods  are  also  implemented  as  numerical  codes.  Since 
a  constant  strength  normal  dipole  is  equivalent  to  a  pair  of  point  vortices  at  the 
panel  edges,  a  distribution  of  point  vortices  is  used  instead  of  the  piecewise  constant 
dipole  distribution.  In  the  mixed  source  and  vortex  method,  point  sources  are  also 
distributed  on  the  same  location  as  the  vortices.  The  collocation  point  is  selected 
on  the  exact  foil  surface  instead  of  the  midpoint  of  each  panel.  An  explicit  Kutta 
condition  is  imposed  at  the  trailing  edge  by  setting  the  tangential  velocity  equal  to 


29 


Total  potential, 


Figure  2.5;  Convergence  characteristics  of  the  total  potential  distribution  calculated 
by  the  total  potential  method.  (NACA  66  +  a=0.8,  t/c=0.04,  f/c=0.02,  a=1.5  deg) 
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Figure  2.6:  Convergence  characteristics  of  the  pressure  distribution  calculated  by  the 
total  potential  method.  (NACA  66  +  a=0.8,  t/c=0.04,  f/c=0.02,  a=1.5  deg) 
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Table  2.1:  Effect  of  the  number  of  chordwise  panels  on  computed  lift  and  drag  co¬ 
efficient  for  the  two-dimensional  foil  section.  (NACA  66  mod.  a=0.8,  t/c=0.04, 
f/c=0.02,  a=1.5  deg.) 


Pert,  potential 

Total  ] 

potential 

N 

Cl 

Cd 

Cl 

(C'L)r 

20 

0.394 

0.0027 

0.389 

0.374 

0.395 

40 

0.410 

0.0019 

0.407 

0.407 

0.413 

80 

0.420 

0.0006 

0.419 

0.401 

0.422 

160 

0.425 

0.0002 

0.425 

zero.  The  direction  of  the  tangential  velocity  is  defined  to  be  normal  to  the  bisector 
of  the  upper  and  lower  surface  at  the  trailing  edge.  With  this  arrangement  of  point 
vortices  and  sources,  both  methods  give  exact  pressure  values  at  the  collocation  points 
for  a  circular  section,  regardless  of  the  number  of  panels.  But  for  a  thin  section  the 
pressure  distributions  by  both  methods  give  erroneous  results. 

Pressure  distributions  by  the  velocity  field  panel  methods  are  shown  in  Figure  2.7 
with  40  panels  for  the  same  foil  geometry.  Both  the  vortex  method  and  the  mixed 
vortex  and  source  method  give  erroneous  pressure  distributions,  even  though  those 
results  converge  to  the  correct  pressure  distribution  as  the  number  of  the  panels  is 
increased.  One  reason  for  this  erroneous  result  is  that  Equation  2.14  is  a  Fredholm 
integral  equation  of  the  first  kind,  and  the  influence  coefficient  matrix  resulting  from 
this  equation  is  not  diagonally  dominant,  which  makes  the  solution  unstable.  Another 
explanation  might  be  the  fact  that  the  kernel  of  the  equation  for  the  velocity  field 
methods  is  one  order  more  singular  than  that  for  the  potential  field  methods. 

From  the  comparison  of  the  results  by  different  panel  methods,  the  following  con¬ 
clusions  are  made: 

•  While  practically  all  methods  work  well  for  thick  sections,  the  potential  mpthod 

is  substantially  more  accurate  for  very  thin  sections.  This  is  particularly  im- 
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Figure  2.7:  Comparison  of  the  pressure  distributions  by  potential  method  and  the 
velocity  methods.  (NACA  66  +  a=0.8,  t/c=0.04,  f/c=0.02,  0=1.5  deg) 
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portant  for  marine  propellers,  where  thickness/chord  ratios  typically  vary  from 
twenty  percent  at  the  root  to  as  little  as  two  percent  near  the  tip. 

•  The  influence  coefficients  for  the  potential  induced  by  unit  source  and  dipole 
distributions  are  one  order  less  singular  than  the  corresponding  influence  coeffi¬ 
cients  for  the  velocity.  As  a  result,  potential  based  methods  are  expected  to  be 
less  sensitive  to  errors  caused  by  irregular  panelling. 

•  The  computation  of  the  panel  influence  coefficients,  which  is  a  major  contributor 
to  the  total  computing  effort,  is  faster  for  a  potential  method  than  for  a  velocity 
method. 

•  Since  the  potential  influence  coefficients  are  scalar  quantities,  the  required  stor¬ 
age  for  the  potential  method  is  one  third  as  great  as  the  storage  for  a  velocity 
method. 

•  Both  the  perturbation  potential  and  the  total  potential  methods  give  accurate 
results  with  similar  convergence  characteristics.  However,  in  the  case  of  non- 
uniform  inflow  velocity,  the  inflow  velocity  potential  is  difficult  to  define,  hence 
making  the  total  potential  method  difficult  to  apply. 

•  For  a  marine  propeller  application,  the  perturbation  potential  method  is  chosen 
and  will  be  implemented  for  three  dimensional  cases  in  the  following  chapters. 

•  In  principle  a  high-order  panel  method  should  be  more  accurated  than  a  low- 
order  method  for  a  given  number  of  elements.  However  a  properly  formulated 
low-order  method  can  achieve  similar  accuracy  by  increasing  the  number  of 
elements.  Since  the  high-order  method  requires  more  computational  cost  per 
element,  the  choice  is  determined  by  the  trade-off  between  the  cost  per  element 
and  the  accuracy  for  a  given  number  of  elements.  Youngren  et.  al.  [24]  claims 
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that  a  low-order  potential  based  panel  method  is  more  efficient  than  a  high-order 
potential  method  in  the  case  of  a  subsonic  flow  and  an  incompressible  flow. 
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Chapter  3 


Numerical  formulation 


By  comparing  the  characteristics  of  various  panel  methods,  the  perturbation  potential 
method  was  chosen  in  Chapter  2  as  the  most  suitable  method  for  the  application  of  the 
marine  propeller  problem.  In  this  chapter,  the  detailed  numerical  implementation  of 
the  perturbation  potential  method  for  three  dimensional  problems  will  be  presented. 

A  discretized  form  of  the  integral  equation  (Equation  2.9)  can  be  applied  to  an 
arbitrary  general  body  in  potential  flow.  The  body  and  wake  surfaces  are  replaced  by 
a  large  number  of  plane  quadrilateral  panels,  and  the  singularity  strength  distribution 
on  these  surfaces  is  approximated  by  a  piecewise  constant  distribution  over  the  panels. 
The  control  point,  where  the  discretized  integral  equation  is  satisfied,  is  selected  as 
the  centroid  of  each  panel.  This  results  in  a  system  of  linear  algebraic  equations  for 
the  unknown  dipole  strengths  which  are  also  potential  values.  From  the  solution  of 
the  linear  system,  any  flow  quantities  of  physical  interest  can  be  calculated. 

Because  our  main  object  of  interest  is  a  marine  propeller  in  steady  flow,  the  nu¬ 
merical  procedure  in  this  chapter  will  be  focused  on  the  case  of  marine  prooeller. 
However  the  numerical  procedure  for  a  general  body  would  be  similar  to  this. 
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3.1  Blade  geometry 


The  propeller  consists  of  K  identical  blades  and  an  axisymmetric  hub.  For  the  steady 
flow  problem,  both  the  geometry  and  the  singularity  distribution  is  repeated  iden¬ 
tically  on  each  blade  and  on  each  inter-blade  segment  of  the  hub.  In  order  to  take 
advantage  of  the  rotational  symmetry  of  the  propeller  problem,  only  ^  portion  of  the 
propeller,  i.e.,  one  blade  and  an  inter-blade  segment  of  the  hub,  is  discretized.  Effect 
of  the  other  portion  of  the  propeller  is  included  in  the  calculation  of  the  influence 
functions,  ais  will  be  explained  in  section  3.5. 

The  propeller  geometry  problem  consists  of  finding  the  Cartesian  coordinates  of 
points  on  the  actual  propeller  surface,  given  the  usual  propeller  geometric  descriptions. 
The  geometry  is  specified  with  respect  to  a  right-handed,  blade-fixed  coordinate  sys¬ 
tem,  with  the  x-cLxis  pointing  downstream  and  the  y-axis  at  some  arbitrary  angular 
orientation  relative  to  the  selected  blade.  The  2-axis  completes  the  right-handed 
system.  Cylindrical  coordinates  {x,r,6)  are  defined  as  usual,  with 

r  =  ^/y2  +  2^  (3.1) 

and  6  being  measured  clockwise  from  the  y-ajcis  when  viewed  looking  downstream. 

The  blade  is  formed  starting  with  a  midchord  line,  which  is  a  space  curve  defined 
parametrically  by  the  radial  distribution  of  skew  and  rake  Xm{r).  By  advancing 

a  distance  along  a  helix  of  pitch  angle  d>(r)  passing  through  the  midchord  line, 

one  obtains  the  blade  leading  and  trailing  edges.  The  nose-tail  line  is  defined  as  a 
helical  line  between  the  leading  and  trailing  edge  at  each  radius.  The  blade  mean 
camber  surface  may  then  be  defined  by  adding  the  camber,  /(r),  at  right  angle  to 
the  nose-tail  line  at  each  radius.  Finally,  thickness  t{r)  is  added  symmetrically  with 
respect  to  the  mean  camber  line  at  each  radius,  again  in  a  cylinder  of  radius  r,  and 
at  right  angles  to  /.  The  coordinate  systems  and  the  notations  described  above  are 
illustrated  in  Figure  3.1. 


37 


Figure  3.1:  Propeller  blade  geometry  notation. 
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Nc  panels  are  distributed  on  each  chordwise  strip  of  the  blade.  The  chordwise  panel 
arrangement  is  the  same  as  that  of  the  two-dimensional  hydrofoil,  and  cosine  spacing 
is  again  adopted.  The  radial  distribution  of  panel  size  is  determined  to  concentrate 
the  elements  toward  the  tip.  In  this  case,  the  radii  for  Mr  panels  are, 

rr„  =  Rf,  +  {R-  Rh)  m  =  1,2,...,M, +  1  (3.2) 

where  R^  is  the  hub  radius,  R  is  the  tip  radius,  and  Mr  is  the  number  of  panels  over 
the  radius.  Panel  arrangement  for  a  three  bladed  propeller  is  illustrated  in  Figures  3.2 
and  3.3,  where  panels  are  distributed  over  the  full  propeller. 

3.2  Hub  geometry 

The  geometry  of  the  hub  is  defined  by  a  profile  curve,  which  can  be  anything  from  a 
constant  diameter  cylinder  to  a  complete  aixisymmetric  body  on  which  the  propeller 
is  mounted.  A  realistic  geometry  may  be  a  semi-infinite  cylinder  on  the  upstream 
side,  with  a  given  fairwater  geometry  downstream  of  the  blades.  The  geometry  of  the 
hub  which  is  used  in  the  sample  calculation  in  Chapter  4  is  shown  in  Figure  3.4.  Here 
the  profile  curve  of  the  fairwater  is  given  as  a  quartic  function  of  the  nondimensional 
distance  from  the  tail  of  the  hub. 

The  geometry  of  the  hub  can  be  determined  by  the  following  input  parameters: 

•  Maximum  hub  radius,  R^. 

•  Computational  ELxial  distance  in  the  upstream  region,  Xu. 

•  Axial  distance  from  the  trailing  edge  of  the  blade  to  the  beginning  of  the  fair- 
water,  Xj. 

•  Axial  length  of  the  fairwater,  Xt. 
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Figure  3.2:  Panel  arrangement  viewed  from  upstream  for  a  three  bladed  propeller 
with  a  hub.  {N,=A0,  A/,  =  10,  Nu  +  Nt,  +  Nci  +  Nt=39,  Me=8) 
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Figure  3.3:  Panel  arrangement  viewed  from  downstream  for  a  three  bladed  propeller 
with  a  hub.  (^^^=40,  Mr=10,  +  Ni,  +  Nd  +  Nt=39,  M$=8) 
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Figure  3.4:  Geomertry  notation  of  the  hub. 
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The  axial  distance  between  the  leading  and  trailing  edge,  Xj  in  Figure  3.4,  is  deter¬ 
mined  from  the  geometry  of  the  blade  at  the  root  section.  Xj  and  Xj  in  Figure  3.4 
are  the  distances  from  the  origin  to  the  upstream  end  and  to  the  downstream  end  of 
the  hub. 

The  number  of  panels  on  the  hub  can  be  determined  by  the  following  integers: 

•  Number  of  circumferential  panels  in  the  non-redundant  region  of  the  hub.  Mg. 


•  Number  of  axial  panels  in  the  upstream  of  the  blades,  N^- 

•  Number  of  axial  panels  between  the  leading  and  trailing  edge  of  the  blade,  Ng. 

•  Number  of  aaial  panels  from  the  trailing  edge  to  the  beginning  of  the  fairwater, 

N,. 

•  Number  of  axial  panels  in  the  fairwater  region,  Nf 


Here,  Ng  is  determined  from  the  blade  panel  arrangement,  Nj  and  Nt  are  determined 
from  the  wake  panel  arrangement,  while  Mg  and  are  given  as  inputs. 

The  arrangement  of  hub  panels  is  selected  to  minimize  the  possible  discretization 
errors  due  to  the  mismatch  of  the  hub  panels  and  the  blade  or  wake  panels.  Over  the 
axial  location  between  the  leading  and  trailing  edge,  the  hub  panels  are  arranged  to 
match  the  blade  panels  at  the  intersections,  and  the  circumferential  panels  along  the 
hub  are  selected  to  have  an  equal  circumferential  angle  of 

The  panelling  on  the  hub  upstream  of  the  propeller  leading  edge  is  purely  helical, 
with  a  pitch  matching  the  root  section  pitch  of  the  propeller.  The  axial  spacing  is 
chosen  to  provide  a  fine  spacing  near  the  propeller  and  a  coarse  spacing  upstream.  In 
this  case,  the  axial  coordinates  of  the  panel  boundaries  are, 


=  -Xj  -h  XuSin( 


7r(n  —  l) 


),  n  —  l,2,...,Xu-t-l. 


(3.3) 
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The  arrangement  on  the  hub  panels  downstream  of  the  trailing  edge  is  similar,  except 
that  the  pitch  is  required  to  match  the  corresponding  pitch  of  the  wake  at  the  hub 
intersection. 

The  axial  coordinates  of  the  intermediate  panels  between  the  blades  must  be  ad¬ 
justed,  particularly  near  the  leading  edge,  in  order  to  avoid  badly  shaped  panels  in 
this  region.  Since  finite  thickness  is  superimposed  to  the  mean  camberline,  the  axial 
coordinate  of  the  panel  vertex  in  the  suction  side  of  the  leading  edge  can  be  less  than 
that  of  the  leading  edge.  If  the  axial  coordinates  of  the  upsteam  panels  were  required 
to  be  the  same  in  the  circumferential  direction,  panels  near  the  leading  edge  could 
actually  turn  inside-out.  To  avoid  this,  the  panel  boundary,  which  intersects  the  suc¬ 
tion  side  of  the  blade  panel  boundary  and  the  upstream  hub  helix  at  the  leading  edge, 
is  selected  to  be  a  bisector  of  the  upstream  helix  and  the  blade  panel  boundary. 

Once  the  vertex  of  the  first  panel  at  the  leading  edge  is  obtained  by  advancing 
the  same  circumferential  interval,  the  axial  coordinates  of  the  intermediate  panels 
along  the  hub  are  obtained  in  a  smooth  manner  as  illustrated  in  Figure  3.5,  where 
the  expanded  plan  view  of  the  hub  panel  arrangement  is  shown. 

The  influence  of  the  upstream  hub  with  semi-infinite  extent  can  be  accounted  for 
if  is  increased  until  no  significant  change  in  the  potential  values  near  the  blades 
is  detected.  However  it  is  observed  that  the  dipole  strength  on  the  far  upstream  hub 
becomes  a  constant  as  the  value  of  is  increased.  While  the  source  panels  in  the  far 
upstream  hub  have  no  influence  because  of  their  zero  strength,  the  dipole  panels  of 
constant  strength  over  the  cylindrical  surface  of  a  semi-infinite  circular  cylinder  can 
be  replaced  by  dipole  panels  over  the  circular  disk  at  the  beginning  of  the  cylinder.  In 
this  way,  the  computational  axial  distance  of  the  upstream  hub,  X^,  can  be  reduced 
significantly. 
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Figure  3.5:  Expanded  plan  view  of  the  hub  panel  arrangement. 
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3.3  Geometry  of  the  trailing  wake 


The  influence  of  the  trailing  wake  on  the  blade  potential  is  a  function  of  both  its 
strength  and  geometry.  The  strength  can  be  related  to  the  radial  distribution  of 
circulation  on  the  blade.  The  geometry,  however,  follows  in  an  indirect  way  from  the 
requirement  that  each  element  of  the  trailing  vorticity  must  be  aligned  with  the  local 
flow.  Since  this  local  flow  depends  on  the  trailing  wake  geometry  in  a  nonlinear  way, 
an  iterative  procedure  must  be  employed.  Greeley  and  Kerwin  '5]  developed  a  wake 
alignment  scheme  which  is  extremely  fast,  yet  capable  of  providing  the  aligned  wake 
geometry.  The  wake  model  in  the  present  work  follows  that  given  in  [5],  and  only  a 
brief  description  will  be  given  in  this  section. 

The  propeller  wake  is  divided  into  two  parts: 

•  A  transition  wake  region  where  the  contraction  and  deformation  of  the  slip¬ 
stream  occurs. 

•  An  ultimate  wake  region  which  is  composed  of  K  concentrated  helical  tip  and 
hub  vortices. 

The  axial  variation  of  the  radii  of  the  trailing  wake  can  be  determined  by  a  limited 
set  of  parameters,  chosen  in  accordance  with  experimental  data. 

•  The  ultimate  radius  of  the  contracted  slipstream,  R^. 

•  The  radius  of  the  hub  vortex  at  the  end  of  the  transition  wake,  Ry,h- 

•  The  length  of  the  transition  wake  region,  Xtw 

•  The  contraction  angle  of  the  tip  vortex  as  it  leaves  the  blade  tip,  6c. 

The  radius  of  the  outermost  trailing  vortex  in  the  transition  wake  is  set  by  a 
smooth  curve  consistent  with  the  above  wake  descriptors. 
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The  radius  of  the  innermost  vortex  is  set  to  that  of  the  fairwater  until  it  becomes 
smaller  than  R^uh-  Near  the  apex  of  the  hub  where  the  radius  of  the  fairwater  is 
smaller  than  the  radius  of  the  innermost  trailing  vortex  is  set  to  R^h-  The  radii 
of  intermediate  trailers  are  obtained  by  an  interpolation  in  the  same  manner  as  was 
described  in  [5].  Figure  3.6  illustrates  the  axial  variation  of  radius  of  a  set  of  trailing 
vortex  elements  derived  by  the  above  description. 

The  pitch  of  the  transition  waJce  must  be  allowed  to  vary  to  match  the  local  fluid 
velocity  calculated  on  the  wake  element.  An  iterative  procedure  is  required  at  this 
stage  since  the  convection  velocities  depend  on  the  geometry  of  the  wake.  In  the 
present  work  a  fully  converged  wake  geometry,  obtained  from  the  lifting  surface  code 
with  the  wake  alignment  scheme  [5],  is  given  as  an  input. 

Since  normal  dipoles  are  employed  to  represent  the  wake,  dipole  panels,  instead 
of  vortex  lines,  are  distributed  on  the  wake  surface  in  the  present  work. 

3.4  Approximation  for  the  induced  potential  of  the 
ultimate  wake 

The  concentrated  tip  and  hub  vortices  in  the  ultimate  wake  region  must,  in  principle, 
extend  downstream  to  infinity.  Because  dipole  panels,  instead  of  vortex  lines,  are  used 
to  represent  the  wake  surface,  helical  strips  of  dipole  panels  also  must  be  extended  to 
infinity. 

Because  the  ultimate  wake  is  located  far  downstream  of  the  blades,  the  induced 
potential  due  to  the  ultimate  wake  panels  can  be  approximated  by  that  of  a  sink  disk 
at  the  beginning  of  the  ultimate  wake  region.  This  will  save  much  computing  effort, 
particularly  for  the  present  method  which  is  based  on  the  potential  field  formulation, 
since  the  induced  potential  due  to  a  dipole  panel  goes  to  zero  at  a  slower  rate  than 
the  induced  velocity  as  the  distance  between  the  field  point  and  the  panel  becomes 
large. 
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Figure  3.6:  Radii  of  the  trailing  vortex  lines.  (A/,  =  10,  i2^=0.83,  X(^  =  1.5, 

<5^  =  15  deg.) 
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Consider  K  equally  spaced  identical  helical  vortex  lines  of  semi-infinite  extent 
with  a  strength  T  and  a  pitch  angle  on  a  circular  cylindrical  surface  of  radius 
as  shown  in  Figure  3.7.  These  vortex  lines  are  connected  to  the  center  line  vortex  via 


radial  vortex  lines  at  the  capping  surface  of  the  cylinder  and  downstream  at  infinity. 
This  system  of  vortex  lines  is  equivalent  to  K  helical  strips  of  normal  dipoles  with 
their  axis  normal  to  the  helical  surface. 

For  a  field  point  far  upstream,  the  influence  (i.e.,  the  induced  velocity  or  potential) 
due  to  the  system  of  these  vortex  lines  can  be  approximated  as  that  of  an  axisym- 
metric  distribution  of  helical  vorticity  on  the  cylindrical  surface  and  a  radial  vorticity 
distribution  on  the  capping  surface  of  the  circular  cylinder,  which  is  connected  to  the 
centerline  vortex.  This  system  of  vorticity  can  be  regarded  as  a  vorticity  representa¬ 
tion  of  a  propeller  with  infinite  number  of  blades.  The  helical  vorticity  has  the  same 
pitch  angle  as  the  vortex  lines  and  has  strength  7  =  • 

This  system  of  the  vorticity  can  be  decomposed  into  two  components  (Figure  3.7): 

1.  A  system  of  vorticity  on  the  cylindrical  surface  with  its  direction  chosen  as 
parallel  to  the  i-axis,  the  radial  vorticity  on  the  capping  surface  of  the  cylinder, 
and  the  centerline  vortex.  The  strength  of  the  helical  vorticity  is  '^s\n4>. 

2.  A  distribution  of  ring  vortices  on  the  circular  cylindrical  surface  with  its  strength 
7  cos  (f>. 

The  system  of  i-vorticity  is  equivalent  to  a  distribution  of  dipoles,  with  their 
axis  in  the  —6  direction,  inside  the  circular  cylinder.  The  induced  potential  at  the 
upstream  field  point  due  to  this  dipole  distribution  is  zero,  because 


(3.4) 


On  the  other  hand,  the  distribution  of  ring  vortices  on  the  circular  cylinder  is 


equivalent  to  a  distribution  of  i-directed  dipoles  inside  the  cylinder.  The  induced 
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Figure  3.7:  Far-field  approximation  of  the  ultimate  wake. 
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potential  due  to  this  dipole  distribution  is  equivalent  to  that  due  to  a  sink  disk  of  the 
same  strength  at  the  capping  surface  of  the  cylinder,  since 

/•JJ.  f2r  roo  Q  _r  /•/?.  r2x  _1  C  t  \ 

i  */o  ''‘‘L  ‘‘^Tx-r^L  ‘‘'i  ,/  r'^^- 

capp^ngsurf  ace 


As  a  result,  the  ultimate  wake  surface  is  replaced  by  a  sink  disk  at  the  beginning  of 
the  ultimate  wake  region. 

A  simple  numerical  verification  of  the  above  approximation  was  performed,  and 
the  results  are  summarized  in  Table  3.1. 


Table  3.1:  Comparison  of  the  induced  potential  due  to  the  helical  dipole  strips  and 
that  due  to  the  sink  disk  on  the  capping  surface. 


Field  point 
{x,y,z) 

Helical  dipole  strip,  x 

Sink  disk,  N, 

20  X  10 

80  X  10 

40  X  10 

40  X  20 

40  X  40 

20 

Analytic 

(-1.0, 0,0) 

3.703 

3.748 

3.739 

3.848 

3.905 

3.849 

3.902 

(-1.0, 1,0) 

2.932 

2.974 

2.966 

3.074 

3.131 

3.144 

3.191 

For  the  purpose  of  numerical  verification,  three  helical  strips  are  equally  spaced 
inside  the  circular  disk  and  the  pitch/diameter  of  the  outside  helix  is  given  as  one. 
Each  helical  dipole  strip  is  represented  by  x  N^v  panels,  where  N,x,  is  the  number 
of  panels  in  one  revolution  of  the  helix,  and  Nrev  is  the  number  of  revolutions  of  the 
helix.  The  sink  disk  is  replaced  by  N,  triangular  panels  at  the  beginning  of  the  helical 
strip.  Axial  distance  between  the  field  points  and  the  beginning  of  ultimate  wake  is 
given  as  1.0. 

Induced  potential  due  to  the  helical  dipole  strips  are  calculated  by  summing  the 
individual  contributions  from  each  panel.  The  analytic  result  of  the  induced  potential 
due  to  the  sink  disk  of  the  strength  where  K  =  3,  T  =  =  1,  is  also 

given  as  a  reference. 
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The  first  field  point  is  located  on  the  extension  of  the  center  line  of  the  helix. 
Since  induced  potential  divided  by  the  number  of  the  helical  strips  is  constant  for 
any  number  of  helical  strips  on  the  field  points  along  the  center  line,  the  replacement 
of  the  ultimate  wake  by  the  sink  disk  does  not  involve  any  approximation.  This  is 
numerically  confirmed  in  Table  3.1.  The  induced  potential  due  to  the  helical  dipole 
strip  on  the  first  field  point  converges  to  the  analytic  induced  potential  due  to  the 
sink  disk  as  and  Nrev  are  increased.  It  is  also  shown  that  the  result  by  20  sink 
panels  on  the  capping  surface  of  the  cylinder  is  as  accurate  as  that  by  40  x  20  dipole 
panels  on  each  helical  strip,  which  is  equal  to  2400  panels. 

The  second  field  point  is  located  at  the  same  radius  as  the  wake  cylinder,  where 
error  due  to  the  farfield  approximation  is  maximum.  The  difference  between  the  result 
using  40  X  40  dipole  panels  and  the  analytic  result  from  the  sink  disk  may  be  regarded 
as  an  error  introduced  by  the  farfield  approximation.  Again,  20  _.nk  panels  on  the 
capping  surface  of  the  cylinder  is  sufficient  to  give  a  very  accurate  result. 

3.5  Discretization  of  the  singularity  distribution. 

The  solution  of  the  boundary  value  problem  consists  of  determining  the  strength  of 
singularities  representing  the  propeller  and  the  trailing  wake,  so  as  to  satisfy  the 
integral  equation  (Equation  2.9).  The  continuous  distribution  of  singularities  on  the 
propeller  and  wake  surface  is  approximated  by  a  stepwise  distribution  over  the  quadri¬ 
lateral  panels.  Since  the  strengths  of  sources  are  prescribed  by  the  kinematic  boundary 
condition  (Equation  2.2),  only  the  dipole  strengths  are  to  be  determined. 

As  described  in  section  3.1,  surface  panels  are  distributed  only  on  the  selected 
blade  and  the  inter-blade  segment  of  hub.  Effects  of  the  other  portions  of  propeller 
are  included  by  summing  the  influence  functions  calculated  by  rotating  the  control 
points  an  angle  ^  about  the  propeller  axis  until  all  blades  are  accounted  for.  This  is 
equivalent  to  having  equal  panelling  on  all  blades. 
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The  four  vertices  of  the  panel  are  located  on  the  exact  surface  of  the  propeller  and 
the  wake,  which  in  general  are  not  on  a  plane.  The  vertices  of  the  plane  panel  are 
constructed  as  follows.  Assuming  that  four  offset  points  are  specified  on  the  exact 
surface,  adjacent  offset  pairs  are  connected  by  straight  segments.  Their  midpoints 
can  be  shown  to  lie  on  a  plane.  The  plane  vertices  are  then  obtained  by  projecting 
the  offset  points  on  that  plane.  This  construction  minimizes  in  a  mean-square  sense 
the  distance  of  the  offsets  from  the  quadrilateral  plane.  In  addition,  the  sides  of 
adjacent  panels  only  contact  each  other  at  a  single  point,  so  their  interior  surfaces  do 
not  intersect. 

The  influence  functions  for  the  quadrilatral  source  and  dipole  panels  are  computed 
using  the  formulation  developed  by  Newman  [21].  The  induced  potential  is  computed 
exactly  for  nearby  panels,  while  it  is  approximated  by  a  multipole  expansion  for  more 
distant  panels.  Finally,  the  panels  which  are  sufficiently  distant  are  treated  as  point 
sources  and  dipoles.  This  is  really  the  heart  of  any  panel  code,  and  it  is  essential  that 
the  computer  code  for  the  influence  function  be  both  robust  and  efficient. 

The  dipole  strength  of  each  panel,  which  is  also  a  potential  value  on  the  panel,  is 
determined  by  satisfying  a  discretized  form  of  the  integral  equation  (Equation  2.9)  on 
the  control  point,  which  is  chosen  as  the  centroid  of  the  panel.  Adopting  piecewise 
constant  singularity  strengths  on  the  plane  quadrilateral  panel,  the  integral  equation 
can  be  expressed  as  a  system  of  linear  equations. 

"g"  +  E  "'.-(A*),  =  i  =  1,2 . (3,6) 

;  =  1  m=l  y=:l 

where  Npanti  is  the  total  number  of  panels  on  the  blade  and  the  hub,  Mr  is  the  number 
of  radial  panels  on  the  blade,  Aj  is  the  induced  potential  at  the  i-th  control  point  due 
to  the  j-th  dipole  panel  of  unit  strength,  S,j  is  the  induced  potential  at  the  i-th  control 
point  due  to  the  j-th.  source  panel  of  unit  strength,  Wi^  is  the  induced  potent-al  at 

the  :-th  control  point  due  to  the  M-th  streamwise  dipole  strip  in  the  wake,  and  (|^)j 
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is  the  source  strength  at  the  j-th  panel  given  by  the  boundary  condition.  Since  the 
strength  of  the  M-th  streamwise  dipole  strip,  can  be  related  to  the  potential 

values  on  the  blade  by  a  numerical  Kutta  condition,  Equation  3.6  is  sufficient  to 
determine  the  unknown  potential  values. 

3.6  Kutta  condition 

A  numerical  Kutta  condition  must  be  imposed  to  specify  the  circulation  around  the 
body,  which  is  equal  to  the  potential  jump  in  the  wake  surface.  An  approximate 
Kutta  condition  for  the  potential  bcised  panel  method,  which  was  first  introduced  by 
Morino  [20],  required  that  the  strength  of  the  dipole  sheet  in  the  wake  be  equal  to 
the  difference  in  the  value  of  the  dipole  strengths  of  the  two  panels  adjacent  to  the 
trailing  edge. 

Morino ’s  Kutta  condition  gives  accurate  results  for  thin  two-dimensional  foil  sec¬ 
tions  with  small  trailing  edge  angles.  However,  it  is  found  that  this  form  of  Kutta 
condition  contains  a  fundamental  error  when  the  free  stream  contains  a  component 
in  the  direction  of  a  line  connecting  the  control  points  of  the  two  trailing  edge  panels. 
Therefore,  for  three  dimensional  problems  with  significant  cross  flows,  a  new  form  of 
Kutta  condition  is  needed  to  satisfy  the  zero  loading  condition  at  the  trailing  edge. 

For  two-dimensional  problems  where  no  cross  flow  component  exists,  a  Kutta 
condition  can  be  the  requirement  that  the  potential  jump  in  the  wake  should  be  equal 
to  the  difference  in  the  total  potential  values,  instead  of  the  perturbation  potential 
values,  of  the  upper  and  lower  panels  at  the  trailing  edge. 

{M),vakc  =  -  d>'  +  f/oo  •  r, (3.7) 

where  f(,j.  is  the  vector  between  the  the  control  points  of  the  two  trailing  edge  panels. 
The  correction  term  to  Morino’s  Kutta  condition,  Uoo  ■  represents  the  potential 
jump  between  trailing  edge  control  points  due  to  the  inflow  velocity  and  vanishes  for 
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cusped  foil  sections  as  the  upper  and  lower  control  points  approach  the  trailing  edge. 

Figures  3.8  and  3.9  show  the  perturbation  potential  and  the  pressure  distribution 
for  a  circular  cylinder  at  90  degree  angle  of  attack,  calculated  by  the  two-dimensional 
panel  code  with  and  without  the  additional  correction  for  the  potential  jump.  Without 
this  correction,  the  circulation  converges  to  about  half  the  correct  value  regardless  of 
the  panel  density.  With  the  correction  term,  the  potential  and  pressure  distributions 
recover  the  analytic  results.  But  it  should  be  noted  that  the  differences  between 
the  results  obtained  by  the  original  Morino  Kutta  condition  are  very  small  for  thin 
two-dimensional  sections. 

For  three  dimensional  problems  w’ith  significant  cross  flows,  this  is  still  not  satis¬ 
factory.  For  the  region  near  the  tip,  where  the  cross  flow  component  around  the  tip 
from  the  pressure  side  to  the  suction  side  of  the  blade  is  prominent,  the  result  by  the 
modified  Morino  Kutta  condition  shows  a  spurious  non-zero  loading  at  the  trailing 
edge.  A  new  form  of  Kutta  condition  is  required  which  enforces  the  zero  loading 
condition,  as  explained  below. 

Figure  3.10  shows  the  flow  at  the  trailing  edge  near  the  tip.  Total  velocity  on 
the  pressiire  (lower)  side,  has  an  outward  component,  while  that  on  the  suction 
(upper)  side,  V'^ ,  has  an  inward  component  due  to  the  cross  flow  around  the  tip.  They 
should  have  the  same  magnitude  to  satisfy  the  zero  loading  condition,  since  the  static 
pressure  in  steady  flow  is  the  same  on  the  upper  and  the  lower  side  of  the  wing  at  the 
trailing  edge. 

Define  a  local  coordinate  system,  where  s  is  chosen  as  the  direction  of  mean  ve¬ 
locities  of  the  upper  and  lower  surfaces,  and  n  is  chosen  normal  to  s  toward  the  tip. 
The  trailing  edge  velocities  can  be  decomposed  into  into  s,  n  directions,  and  for  a 
zero  loading  at  the  trailing  edge,  the  following  equations  must  be  satisfied. 


=  vl 

(3.8) 

1 

II 

(3.9) 
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Figure  3.8:  Perturbation  potential  distribution  for  a  circular  cylinder  at  90  degree 
angle  of  attack. 
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:  Pressure  distribution  for  a  circular  cylinder  at  90  degree  angle  of  attack. 
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Figure  3.10:  Magnified  flow  at  the  trailing  edge  near  the  tip. 

where  and  are  s,n  components  of  and  respectively. 

In  general,  the  calculated  velocities  at  the  trailing  edge  using  the  modified  Morino 
Kutta  condition  do  not  have  the  same  magnitude.  As  a  result,  the  pressure  values 
on  the  upper  and  lower  surfaces  do  not  match.  This  non-zero  loading  at  the  trailing 
edge  becomes  large  toward  the  tip,  where  the  three-dimensional  cross  flow  component 
becomes  large. 

In  the  present  work,  an  explicit  pressure  Kutta  condition  is  employed,  which  re¬ 
quires  the  pressures  on  the  last  panels  at  the  trailing  edge  be  equal.  The  potential 
jump  in  the  wake  surface  is  expressed  as, 

(A(^)u;a>:c  =  ~  +  ^oo  *  U.t.  +  (3.10) 

where  (AC7p)^  ^  is  the  difference  of  the  pressure  values  of  the  upper  and  lower  pan-  ls  at 

trailing  edge  and  /T  is  a  parameter  to  be  found  to  ensure  the  zero  loading  condition. 
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Due  to  the  nonlinear  aspect  of  the  pressure  Kutta  condition,  an  iterative  solution 
algorithm  is  employed. 

In  the  numerical  implementation,  the  potential  jump  at  the  m-th  dipole  wake  strip 
is 

+  K{AC,Y^-'\  (3.11) 


where  and  (f>]^  are  the  potential  values  of  the  upper  and  lower  panels  at  the  trailing 
edge.  Then  Equation  3.6  becomes 


panel 


Mr 


^panel 


d(j). 


j=i 


m=l 


;=i 


dn 


Mr 


-  J2^irn[{Uoo-r,.,.)mK{AC,)l^-%  i  =  1,2,  .  .  .  ,  N.^.U) 


m=l 


During  the  iteration  process,  the  value  of  ACp  is  obtained  from  the  previous  iteration, 
and  the  value  of  K  is  determined  by  the  Newton-Raphson  method.  However,  the 
influence  coefficient  matrix  is  unchanged  in  the  process.  The  initial  guess,  (A(;5)^^,  is 
taken  by  setting  the  value  of  K  to  zero. 

On  the  other  hand,  the  difference  of  the  n-direction  velocities  at  the  trailing  edge 
panels  in  Figure  3.10  is  related  to  the  magnitude  of  the  trailing  vorticity.  The  trailing 
vorticity  is  in  the  s-direction  and  the  magnitude  of  it  must  be  equal  to  the  spanwise 
derivative  of  the  circulation; 


■-r=RI  =  lV“-V'.:|cos,i  =  -— ,  (3.13) 

where  t/;  is  the  angle  between  n-  and  77-coordinates.  This  is  the  compatibility  condition, 
which  is  a  measure  of  the  consistency  of  the  panel  method. 


3.7  Linear  system  solution 


The  unknown  dipole  strengths  are  determined  by  the  solution  of  the  system  of  linear 

equations.  For  a  small  number  of  unknowns,  such  as  two-dimensional  problems  or 

59 


simple  three-dimensional  problems,  Gauss  reduction  can  be  used  at  the  outset  to 
decompose  the  matrix  into  lower  and  upper  triangular  forms  for  subsequent  solutions 
with  different  right  hand  sides.  But  for  complex  three-dimensional  geometry,  the  time 
required  for  the  Gauss  reduction  method  would  become  prohibitive. 

Since  several  iterations  are  required  to  satisfy  the  pressure  Kutta  condition,  an  ef¬ 
ficient  matrix  solver  is  an  essential  feature  of  the  present  panel  method.  In  the  present 
work,  an  accelerated  iterative  matrix  solver  developed  by  Clark  [4]  is  employed.  This 
method  is  found  to  converge  very  rapidly  for  the  kind  of  influence  coefficient  matrix 
encountered  with  the  present  potential  method,  and  requires  very  small  computing 
time  compared  to  the  Gauss  reduction  method  even  for  a  linear  system  with  small 
number  of  unknowns. 

As  the  number  of  unknowns  becomes  large,  the  computer  memory  to  store  the  full 
influence  coefficient  matrix  increases  in  proportion  to  where  N  is  the  number  of 
unknowns.  Instead  of  storing  all  of  the  matrix  elements  into  the  computer  memory, 
they  are  saved  in  the  outside  storage  and  retrieved  when  neccesary.  In  this  way,  the 
computer  memory  required  can  be  significantly  reduced  even  though  more  computer 
time  is  required  for  storing  and  retrieving  the  matrix  elements. 

Table  3.2  provides  a  comparison  of  the  computing  times  of  the  Gauss  reduction 
method  and  the  accelerated  iterative  method  using  a  DEC  Microvaix  2.  Here,  the 
computing  times  required  for  memory  resident  Gauss  reduction  are  given  for  the  424 
and  724  panel  cases,  while  those  for  840  and  1680  panels  include  the  time  required 
for  external  storage  and  retrieval  of  the  matrix  elements.  The  computing  time  for  the 
matrix  solver  is  for  one  iteration  of  the  pressure  Kutta  condition.  The  effectiveness  of 
the  iterative  solver  is  clearly  shown  here.  The  accelerated  matrix  solver  will  be  used 
for  all  the  sample  calculations  in  Chapter  4. 

Table  3.2  also  tabulates  the  required  computing  times  for  the  calculation  of  the 
influence  coefficient  matrix.  Here  again  the  computing  times  for  840  and  1680  panel 
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cases  include  the  time  required  to  store  the  matrix  elements. 


Table  3.2:  Comparison  of  the  computing  times  of  the  influence  function  calculation 
and  the  matrix  solver. 


N 

Inf.  coefficient 

Matrix  solution 

Gauss 

Iterative 

424 

330  sec 

330  sec 

30  sec 

704 

800 

1510 

77  1 

840 

2100 

180 

1680 

5300 

670 

3.8  Calculation  of  velocities,  pressures,  forces  and 
moments. 

Once  the  potential  values  have  been  determined,  surface  velocities  can  be  calculated 
either  by  numerical  differentiation  of  the  potential  or  by  direct  calculation  of  the 
source  and  dipole  panel  influence  functions.  The  latter  approach  is  found  to  be  not  as 
successful  since  the  velocity  influence  functions  are  more  singular  and  therefore  more 
sensitive  to  the  position  of  the  control  point  within  each  panel. 

A  local  second  order  distribution  of  the  perturbation  potential  is  assumed  on  the 
five  panel  centers  (i.e.,  a  central  panel  and  its  immediate  neighbors),  and  a  local 
tangential  perturbation  velocity  is  obtained  by  differentiation.  The  total  tangential 
velocity  is  obtained  by  a  vector  sum  of  the  perturbation  velocity  and  the  undisturbed 
inflow  velocity. 

As  shown  in  Figure  3.11,  where  u,  u-coordinates  are  formed  by  connecting  the 
midpoints  of  the  sides  of  a  panel  and  f ,  r^-coordinates  are  chosen  as  local  orthogonal 
coordinates,  the  calculated  perturbation  velocites  in  the  u,u  direction  are  not  at  right 

angles  to  each  other.  Each  velocity  component  is  a  projection  of  the  velocty  vector 
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Figure  3.11:  Summation  of  the  velocity  component. 


onto  the  u,v  directions,  respectively.  The  ry-velocity  component  is  calculated  as 

d<f>  sin  t/) 

^  1  (3.14) 

OT]  cos  rp  ' 

where  xp  is  the  angle  between  the  r?  and  v-coordinate.  The  perturbation  velocity  and 
the  inflow  velocity  are  added  in  the  local  orthogonal  coordinates. 

The  pressure  coefficient  is  calculated  using  Bernoulli’s  equation,  and  the  non- 
dimensional  pressure  coefficient  is  deffned  as, 


Cp  = 


P-  Poc 

\pUl 


(3.15) 


Finally,  total  forces  and  moments  are  obtained  by  summation  of  individual  panel 
force  vectors.  In  order  to  obtain  practically  useful  results  which  can  be  compared  with 
experiments,  a  viscous  drag  correction  is  needed.  On  each  panel,  a  viscous  friction 

force  is  added  in  the  direction  of  the  inflow  velocity.  The  viscous  friction  coefficient 

62 


is  given  as  0.003  for  the  calculation  at  the  design  condition  and  that  at  the  off-design 
condition  is  determined  from  the  two-dimensional  section  drag  coefficient  curves. 
Integrated  thrust  and  torque  coefficients  are  calculated  by 

^  pn'^D*  ^  (316) 


(3.17) 


where  x,-,  t/j,  2,-,  (nj),-,  (n^),-,  (n*),-  and  A,-  are,  respectively,  control  point  coordinates, 
unit  normal  vector  components,  and  area  for  i-th  panel.  These  sums  can  be  performed 
separately  over  blades  and  hub.  The  propeller  efficiency  is  defined  as 


27r  Kq  ' 


(3.18) 
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Chapter  4 

Numerical  results 


Three  separate  codes  for  different  problems,  i.e.  a  wing,  an  axisymmetric  duct  and  a 
propeller  problem,  are  implemented  in  order  to  verify  the  theory  given  in  the  previous 
chapters. 

Instead  of  writing  a  general  purpose  code  which  would  take  coordinates  of  the 
panel  vertices  as  inputs,  each  code  is  written  for  a  specific  geometry  with  specific 
inputs  which  are  defined  in  a  natural  way  for  the  geometry.  In  this  way  each  code  is 
kept  to  a  reasonable  size,  and  the  possible  occurence  of  fatal  errors  due  to  inconsistent 
inputs  can  be  minimized.  A  perspective  view  of  the  configuration  with  the  hidden 
lines  removed  can  also  be  plotted  in  order  to  check  the  geometry  before  calculating 
the  influence  coefficients. 

Numerical  results  by  each  code  will  be  presented  in  separate  parts  of  this  chapter. 
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Part  I 

Wing  problem 
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The  wing-body  configuration  is  regarded  as  a  two  bladed  propeller  with  infinite 
pitch,  and  accordingly  the  geometric  part  of  the  code  for  the  wing  problem  is  written 
as  a  special  case  for  a  propeller  problem.  An  option  for  a  fuselage  is  given  in  order  to 
calculate  the  performance  of  the  wing-body  or  wing-only  configuration. 

However,  the  wing  problem  possesses  a  lateral  symmetry  with  respect  to  the  center 
plane,  while  the  propeller  problem  is  rotationally  symmetric  with  a  separating  angle 
of  Reduction  of  the  system  of  equations  is  obtained  by  exploiting  the  symmetry. 
A  frozen  wake  geometry  is  assumed  on  the  surface  formed  by  the  extensions  of  the 
nose-tail  lines  of  the  wing  sections. 

4.1  Ellipsoid  at  zero  angle  of  attack 

The  first  example  is  chosen  to  verify  the  code  by  comparing  the  results  with  the  known 
analytic  results  for  a  simple  geometry.  The  example  is  an  ellipsoid  at  zero  angle  of 
attack  with  semi-lengths  of  its  axes  given  as  a=b=l  and  c=0.1.  The  coordinate  system 
and  the  panel  arrangement  with  20  chordwise  and  10  spanwise  panels  are  shown  in 
Figure  4.1,  where  cosine  spacing  is  chosen  in  both  directions  to  account  for  the  rapid 
changes  in  the  geometry  and  the  singularity  strength  near  the  leading  and  trailing 
edges  and  the  tip. 

The  analytic  expression  of  the  perturbation  velocity  potential  on  the  surface  of 
the  ellipsoid  is  given  in  Lamb  [16],  and  that  of  the  surface  velocities  in  the  desired 
directions  (u,v-directions  in  Figure  4.1)  is  given  in  Appendix  C. 

The  chordwise  potential  and  the  u-velocity  are  compared  with  the  analytic  solu¬ 
tions  in  Figures  4.2  and  4.3  at  four  spanwise  locations.  The  abscissa  is  taken  as  the 
nondimensional  chordwise  parameter  s,  which  is  zero  at  the  lower  trailing  edge,  0.5 
at  the  leading  edge,  and  1.0  at  the  upper  trailing  edge.  The  agreement  is  excellent 
except  for  the  velocities  at  the  leading  and  trailing  edges. 

Since  the  velocity  is  calculated  by  a  difference  formula  from  the  velocity  poten- 
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Figure  4.1:  Coordinate  system  and  panel  arrangement  of  the  ellipsoid. 


67 


tial,  a  small  error  in  the  potential  distribution  would  be  magnified  in  the  velocity 
distribution.  These  magnified  errors  in  the  velocity  distribution  at  y/b=0.980  are 
shown  in  Figure  4.3,  but  the  effect  of  these  errors  is  localized.  The  discrepancy  at  the 
leading  and  trailing  edges  is  due  to  the  error  in  calculating  the  distance  between  the 
control  points,  which  increases  near  the  rounded  edges.  This  error  would  decrease  at 
the  trailing  edge  if  a  conventional  section  were  chosen  instead  of  the  elliptic  section. 
The  velocity  at  the  Icist  panel  at  the  trailing  edge  is  linearly  extrapolated  from  the 
previous  two  values. 

The  spanwise  distribution  of  the  potential  and  the  v-velocity  are  illustrated  in 
Figures  4.4  and  4.5  at  three  chordwise  locations.  The  pressure  distribution  is  cal¬ 
culated  by  Bernoulli’s  equation  using  the  (u,v)  velocities,  as  described  in  Chapter  3. 
Figures  4.6  and  4.7  provide  a  comparison  of  the  chordwise  pressure  distributions  by 
the  present  method  with  the  analytic  solution  at  four  spanwise  locations.  Again,  the 
agreement  is  excellent  except  for  those  very  close  to  the  tip. 

4.2  Circular  wing  at  finite  angle  of  attack 

The  next  example  is  a  circular  planform  wing  with  an  NACA  four  digit  section  at 
an  angle  of  attack  of  0.1  radian.  Thickness  effect  for  the  wing  is  investigated  for  the 
cases  of  two  different  thickness/chord  ratios  of  0.01  and  0.05.  The  analytic  result  for 
the  lifting  circular  wing  with  zero  thickness  was  given  by  Jordan  (12|,  which  can  be  a 
standard  of  comparison  for  the  present  method  as  the  thickness  goes  to  zero. 

The  circular  planform  wing  is  particularly  important  for  the  propeller  problem, 
because  a  magnified  flow  near  the  propeller  tip  is  similar  to  the  flow  around  the 
circular  wing. 

The  coordinate  system  and  panel  arrangement  are  the  same  as  those  of  the  ellipsoid 
case  in  Figure  4.1,  only  the  section  shape  is  the  NACA  four  digit  series  instead  cf  the 

elliptic  section.  Both  the  chordwise  and  spanwise  panelling  are  cosine  spacing. 
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Figure  4.2:  Chordwise  potential  distribution  of  the  ellipsoid.  {a=b=l,  c=0.1) 


69 


a.  la 


T 


T 


r 


T 


-^3  0 . 0O 

<01*0 


r 


Ellipsoid 

(a=b=l,c=0.1,Ar,=40,M,=20) 


—  analytic 
+  r/R=0.039 
0  r/R=0.980 


r 


0.00 


0.00 


0.40  0.60 

s 


0.60 


Figure  4.3;  Chordwise  distribution  of  the  u-velocity  component  of  the  ellipsoid. 
(a=b=l,  c=0.l) 
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Figure  4.5:  Spanwise  distribution  of  the  v-velocity  component  of  the  ellipsoid 
(a=b=l,  c=0.1) 
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Figure  4.6:  Chordwise  pressure  distribution  of  the  ellipsoid.  (a=b=l,  c=0.1) 
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A  series  of  computations  is  performed  for  the  circular  wing  with  an  NACA  0001 
section.  The  radial  circulation  distributions  for  varying  numbers  of  chordwise  panels 
(Nc)  are  presented  in  Figure  4.8  with  Jordan’s  result  for  zero  thickness  circular  wing. 
The  effect  of  the  number  of  spanwise  panels  {Mr)  on  the  circulation  is  shown  in 
Figure  4.9.  For  the  circular  wing  case,  the  spanwise  convergence  is  relatively  slow 
compared  to  the  chordwise  convergence,  which  is  a  typical  characteristic  of  low  aspect 
ratio  wings.  The  results  with  40  chordwise  and  40  spanwise  panels  are  considered 
converged,  and  the  calculation  is  performed  with  this  number  of  panels  hereafter. 

As  described  in  Section  3.6,  the  pressure  Kutta  condition  is  applied  for  this  lifting 
problem.  Figure  4.10  provides  a  comparison  of  the  pressure  distributions  calculated 
by  Morino’s  Kutta  condition  and  by  the  pressure  Kutta  condition.  The  pressure 
differences  between  the  upper  and  lower  surfaces  for  the  circular  wing  with  an  NACA 
0001  section  are  shown  along  the  chordwise  panels  at  a  radial  position  of  r/R=0.916. 
The  pressure  difference  at  the  trailing  edge  by  Morino’s  Kutta  condition  has  a  positive 
value,  which  should  be  set  to  zero  by  the  pressure  Kutta  condition. 

The  pressure  difference  distributions  for  the  wing  with  an  NACA  0005  section 
at  the  same  radial  location  are  shown  in  Figure  4.11.  The  positive  loading  at  the 
trailing  edge  by  Morino’s  Kutta  condition  is  seen  to  be  increased  compared  to  that 
for  the  thinner  section.  Again,  the  positive  loading  is  set  to  zero  by  the  pressure 
Kutta  condition. 

The  non-zero  loading  at  the  trailing  edge  near  the  tip  has  annoyed  many  re¬ 
searchers  in  this  field  [15], [24].  When  the  usual  Morino’s  Kutta  condition  is  applied,  a 
positive  loading  at  the  trailing  edge  near  the  tip  results  for  a  wing  with  large  negative 
sweep  angles,  while  a  negative  loading  results  for  a  wing  with  large  positive  sweep 
angles.  This  can  be  explained  as  follows. 

Figure  4.12  shows  a  magnified  flow  at  the  trailing  edge  near  the  tip  for  a  wing  with 
large  negative  sweep  angles.  A  tip  cross  section  profile  parallel  to  the  trailing  edge 
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Figure  4.8:  Effect  of  the  chordwise  number  of  panels  on  the  circulation  distribution 
of  the  circular  wing.  (t/c=0.01,  a=5.73  deg.) 
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Figure  4.9:  Effect  of  the  spanwise  number  of  panels  on  the  circulation  distribution  of 
the  circular  wing.  (t/c=0.01,  a=5.73  deg.) 
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Figure  4.10:  Effect  of  the  pressure  Kutta  condition  on  the  chordwise  distribution  of 
the  pressure  difference  of  the  circular  wing.  (t/c=0.01,  a=5.73  deg.,  r/R=0.916) 
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Figure  4.11:  Effect  of  the  pressure  Kutta  condition  on  the  chordwise  distribution  of 
the  pressure  difference  of  the  circular  wing.  (t/c=0.05,  a=5.73  deg.,  r/R=0.916) 
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Figure  4.12:  Magnified  flow  at  the  trailing  edge  near  the  tip. 

is  also  shown.  The  inflow  velocity  can  be  decomposed  into  the  x-  and  2-components. 
Due  to  the  2-component  of  inflow  velocity,  a  cross  flow  around  the  tip  results,  which 
is  locally  similar  to  the  flow  around  a  wedge.  The  cross  flow  component  is  outboard 
on  the  lower  surface  and  inboard  on  the  upper  surface.  The  magnitude  of  the  cross 
flow  velocity  becomes  larger  for  a  thicker  section. 

The  cross  flow  is  added  to  the  i-component  inflow  velocity.  With  a  negative  sweep 
angle,  the  magnitude  of  the  resulting  velocity  on  the  upper  surface  is  larger  than  that 
on  the  lower  surface,  which  results  in  positive  loading  at  the  trailing  edge  near  the 
tip.  Since  the  trailing  edge  region  of  the  circular  wing  toward  the  tip  is  geometrically 
similar  to  a  wing  with  negative  sweep  angles,  the  loading  distribution  there  behaves 
similar  to  the  negative  sweep  case.  For  the  case  of  a  wing  with  large  positive  sweep 
angles,  the  situation  is  reversed  and  negative  loading  at  the  trailing  edge  results. 

If  a  wake  alignment  scheme  were  adopted,  the  trailing  vortex  wake  would  leave  the 
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trailing  edge  in  the  mean  direction  of  velocities  on  the  upper  and  lower  surfaces,  and 
the  trailing  wake  would  contract  for  the  negative  sweep  case.  Since  the  contraction 
of  the  wake  would  increase  the  downwash  velocities  at  the  control  points  in  the  inner 
radii  of  the  trailing  vortex,  a  fully  aligned  wake  model  would  alleviate  the  positive 
loading  at  the  trailing  edge.  The  frozen  wake  model  in  the  present  method  does  not 
account  for  the  wake  contraction.  Only  the  positive  loading  should  be  set  to  zero  by 
the  pressme  Kutta  condition. 

As  the  thickness  of  the  circular  wing  increases,  the  magnitude  of  the  cross  flow  com¬ 
ponent  increases  and  the  positive  loading  at  the  trailing  edge  with  Morino’s  Kutta  con¬ 
dition  increases.  In  order  to  satisfy  the  pressure  Kutta  condition,  the  dipole  strength 
in  the  wake  is  reduced  during  the  iteration  process.  As  a  result,  circulation  is  re¬ 
duced  as  the  thickness  is  increased,  which  contradicts  the  trend  for  two-dimensional 
cases.  This  trend  is  shown  in  Figure  4.13,  where  the  radial  circulation  distributions 
are  illustrated  for  the  wing  with  NACA  0001  and  NACA  0005  sections. 

Figures  4.14  through  4.17  show  the  chordwise  distributions  of  the  pressure  differ¬ 
ence  at  four  different  spanwise  locations.  Here  the  results  for  thickness/chord  ratios 
of  0.01  and  0.05  are  plotted  with  Jordan’s  result  for  the  zero  thickness  circular  wing. 
As  one  can  see,  the  result  of  the  0.01  thickness-chord  ratio  wing  agrees  very  well  with 
Jordan’s  result.  This  indicates  the  robustness  of  the  present  method  even  for  the  very 
thin  wing.  Positive  loading  is  shown  near  the  trailing  edge  toward  the  tip  for  the  wing 
with  the  NACA  0005  section.  This  would  decrease  if  a  he  contraction  model  for  the 
trailing  wake  were  adopted. 

4.3  Rectangular  wings  with  different  sweep  angles 

The  next  example  is  chosen  to  study  the  effect  of  sweep  angles  on  the  performance  of 
rectangular  wings.  A  series  of  computations  is  made  for  wings  with  5.9  aspect  ratio, 
8  degree  incidence,  and  varying  sweep  angles  of  -45,  0  and  +45  degrees.  For  each 
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Figure  4.13:  Thickness  effect  on  the  radial  circulation  distribution  of  the  circular 
wing.  {0=5.73  deg.) 
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Figure  4.14:  Chordwise  distribution  of  the  pressure  difference  of  the  circular  wing  at 
r/R=0.098. 
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Figure  4.15t  Chordwise  distribution  of  the  pressure  difference  of  the  circular  wing  at 
r/R=0.693. 
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Figure  4.16:  Chordwise  distribution  of  the  pressure  difference  of  the  circular  wing  at 
r/R=0.916. 
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Figure  4.17:  Chordwise  distribution  of  the  pressure  difference  of  the  circular  wine  at 
r/R=0.977. 
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wing,  the  thickness  effect  is  illustrated  with  the  results  of  the  wings  with  NACA  0001, 
NACA  0005  and  NACA  0012  sections. 

The  panel  arrangements  for  the  wings  with  three  different  sweep  angles  are  shown 
in  Figure  4.18.  Here  again,  cosine  spacing  is  adopted  to  account  for  the  rapid  change 
of  the  singularity  strength  near  the  leading  and  trailing  edges  and  the  tip. 

A  convergence  test  is  performed  for  the  unswept  wing  with  an  NACA  0012  section. 
The  influence  of  number  of  panels  on  circulation  is  illustrated  in  Figures  4.19  and 
4.20.  The  solid  line  in  both  figures  is  the  result  with  80  chordwise  and  20  spanwise 
panels,  which  is  considered  as  converged.  The  chordwise  convergence  is  shown  to  be 
slower  than  the  spanwise  convergence  for  this  relatively  high  aspect  ratio  wing.  For 
a  low  aspect  ratio  wing,  the  convergence  characteristic  would  be  opposite,  cls  shown 
in  the  previous  example  of  the  circular  wing.  Since  the  result  with  40  chordwise  and 
10  spanwise  panels  is  reasonably  accurate  and  requires  much  less  computing  time, 
computation  is  made  with  this  number  of  panels  unless  otherwise  mentioned. 

The  force  calculation  can  be  done  either  by  integration  of  the  element  pressure 
on  each  panel,  or  by  considering  the  energy  far  downstream,  which  is  the  so-called 
Trefftz  plane  wake  integration  method.  The  lift  and  drag  coefficients  calculated  by 
both  methods  are  shown  in  Table  4.1,  for  the  rectangular  wing  with  aspect  ratio 
5.9,  with  an  NACA  0012  section  at  8  degree  incidence.  Nc  represents  the  number 
of  panels  along  the  chordwise  strip  in  the  wing,  while  Mr  represents  the  number  of 
radial  strips.  The  forces  calculated  by  both  methods  converge  to  common  values  as 
the  number  of  panels  increases.  Again  the  result  using  40  chordwise  and  10  spanwise 
panels  is  reasonably  accurate.  The  pressure  integration  method  is  preferred  because 
of  its  direct  physical  interpretation. 

The  thickness  effect  for  the  unswept  wing  is  illustrated  in  Figure  4.21.  For  this 
unswept  wing,  the  circulation  is  increased  as  the  thickness  is  increased,  which  is  the 
same  trend  as  in  the  two-dimensional  Ccise. 
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Figure  4.18:  Panel  arrangements  of  the  rectangular  wings  with  -45,  0,  and  -1-45  degree 
sweep  angles.  (Aspect  ratio=5.9,t/c=0.01,  0.05  and  0.12,  q=8  deg.) 
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Figure  4.19:  Effect  of  the  chordwise  number  of  panels  on  the  circulation  distribution 
of  the  unswept  rectangular  wing.  (Aspect  ratio=5.9,  t/c=0.12,  a=8  deg.) 
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Figure  4.20:  Effect  of  the  spanwise  number  of  panels  on  the  circulation  distribution 
of  the  unswept  rectangular  wing.  (Aspect  ratio=5.9,  t/c=0.12,  a=8  deg.) 
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Figure  4.21:  Effect  of  the  thickness  on  the  spanwise  circulation  of  the  unswept  wing. 
(Aspect  ratio=5.9,  a=8  deg.) 
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Table  4.1:  Effect  of  the  number  of  panels  on  the  lift  and  drag  coefficients  of  the 
unswept  rectangular  wing.  (Aspect  ratio=5.9,  t/c=0.12,  q=8  deg.) 


Nc 

Mr 

(^L)prcjj 

{^L)trefftz 

{(^D)pre$s 

(C'r>)«re//tx 

20 

10 

0.426 

0.410 

0.0189 

0.0145 

40 

10 

0.429 

0.426 

0.0165 

0.0154 

80 

10 

0.432 

0.434 

0.0161 

0.0160 

40 

20 

0.424 

0.422 

0.0160 

0.0148 

80 

20 

0.427 

0.429 

0.0154 

0.0152 

The  pressure  difference  between  the  upper  and  lower  surfaces  is  shown  in  Fig¬ 
ures  4.22  and  4.23  at  three  different  spanwise  locations,  for  thickness/chord  ratios 
of  0.05  and  0.12  respectively.  The  thickness  effect  is  apparent  near  the  leading  edge 
where  a  finite  pressure  peak  is  shown.  The  ability  to  calculate  the  pressure  peak  is 
one  of  the  major  advantages  of  the  present  method  over  the  lifting  surface  theory. 

As  explained  in  Section  3.6,  not  only  the  pressure  Kutta  condition  but  also  the 
compatibility  condition  should  be  satisfied  at  the  trailing  edge.  Hirschel  et.  al.  [ll] 
claimed  that  a  low  order  panel  method,  such  as  the  present  method,  does  not  satisfy 
the  compatibility  condition  near  the  tip.  The  radial  vorticity  distributions  calculated 
by  the  difference  of  the  cross  flow  velocities  between  the  upper  and  lower  surfaces 
are  plotted  with  those  calculated  by  the  spline  differentiation  of  the  circulation  in 
figures  4.24  and  4.25,  for  the  unswept  wings  with  NACA  0005  and  NACA  0012 
sections.  The  agreement  is  clearly  shown,  and  the  compatibility  is  satisfied.  The 
satisfaction  of  the  compatibility  condition  shows  the  self  consistency  of  the  present 
method. 

The  effect  of  thickness  for  the  rectangular  wing  with  -45  degree  sweep  angle  is 
illustrated  in  Figure  4.26,  where  the  distributions  of  the  radial  circulation  are  given  for 
the  varying  thickness/chord  ratios  of  0.01,  0.05  and  0.12.  Since  the  induced  dowr'vash 
due  to  the  trailing  wake  is  large  toward  the  tip  compared  to  the  rectangular  wing,  the 
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Figure  4.22:  Chordwise  distribution  of  the  pressure  difference  of  the  unswept  rectan¬ 
gular  wing.  (Aspect  ratio=5.9,  t/c=0.05,  a=8  deg.) 
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Figure  4.23:  Chordwise  distribution  of  the  pressure  difference  of  the  unswept  rectan¬ 
gular  wing.  (Aspect  ratio=5.9,  t/c=0.12,  a=8  deg.) 
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Figure  4.24:  Spanwise  vorticity  distribution  of  the  unswept  rectangular  wing.  (Aspect 
ratio=5.9,  t/c=0.05,  a=8  deg.) 
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overall  circulation  near  the  tip  is  reduced.  The  effect  of  thickness  very  close  to  the  tip 
is  shown  to  be  the  same  as  that  of  the  circular  wing,  where  the  loading' is  decreased 
as  the  thickness  is  increased.  This  effect  was  explaind  in  the  section  for  the  circular 
wing.  However,  for  the  inner  part  of  the  wing  where  the  cross  flow  component  is 
small,  the  circulation  is  increased  as  the  thickness  is  increased. 

Figure  4.27  shows  the  radial  circulation  distributions  of  the  +45  degree  swept  wing 
with  the  varying  thickness/chord  ratios  of  0.01,  0.05  and  0.12.  Here,  the  downwash  is 
bigger  near  the  inner  part  of  the  wing  due  to  the  positive  sweep,  and  as  a  result  the 
loading  is  reduced.  The  thickness  effect  is  shown  to  have  a  similar  trend  to  that  of 
the  unswept  wing. 

Figure  4.28  provides  a  comparison  of  the  the  local  section  lift  coefficients  by  the 
present  method  with  those  by  other  production  panel  codes,  for  the  unswept  rectan¬ 
gular  wing  with  aspect  ratio  5.9  and  NACA  0012  sections  at  8  degree  incidence.  The 
results  by  the  production  panel  codes  are  taken  from  Margason  et.  al.  [17].  All  of  the 
panel  methods,  including  the  present  one,  over-predict  the  experimental  data.  This  is 
attributable  to  neglecting  the  viscous  effect.  The  result  by  the  present  method  agrees 
well  with  those  by  the  production  codes. 

The  result  for  the  +20  degree  swept  wing  with  aspect  ratio  5.9  and  NACA  0012 
sections  at  8  degree  incidence  is  compared  with  those  by  the  production  codes  in 
Figure  4.29.  Again,  the  result  by  the  present  method  agrees  well  with  those  by  the 
other  methods. 

4.4  Wing-body  configuration 

A  wing-body  configuration  is  examined  to  determine  whether  the  present  method 
can  determine  the  interference  effect  of  a  wing  on  a  fuselage  pressure  distribution. 
Figure  4.30  illustrates  the  wing-body  geometry  selected  for  this  example.  An  exper¬ 
imental  surface  pressure  distribution,  which  was  made  in  the  RAE  8  ft  x  6  ft  wind 
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Figure  4.27:  Radial  circulation  distribution  for  the  +45  degree  swept  wing.  (Aspect 
ratio=5.9,  a=8  deg.) 


99 


LO 


8 


.2 

0 


.1  .2  .3  .4  .5  .6  .7  .8  .9  LO 


r/R 


experiment 
vortex  lattice 
HESS 
VSAERO 
QUADPAN 
•PAN  AIR 

mcaero 

present  metho< 


Figure  4.28:  Spanwise  local  lift  coefficient  distribution  for  the  unswept  rectangular 
wing.  (Aspect  ratio=5.9,  t/c=0.12,  a=8  deg.) 
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Figure  4.29:  Spanwise  local  lift  coefficient  distribution  for  the  +20  degree  swept  rect¬ 
angular  wing.  (Aspect  ratio=5.9,  t/c=0.12,  a=8  deg.) 
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tunnel,  was  taken  from  Treadgold  et.  al.  [22]. 

The  geometry  is  symmetric  with  respect  to  both  the  wing  chord  plane  and  the  zero 
butt  line  plane.  The  cross  section  of  the  wing  is  an  uncambered  RAE  101  section  with 
thickness/chord  ratio  of  0.09.  The  planform  has  an  aspect  ratio  of  6  and  a  mid-chord 
sweep  of  30  degrees.  The  axisymmetric  body  has  a  nose  profile  given  as  a  quartic 
function  of  the  axial  distance  from  the  leading  edge  of  the  body,  which  is  connected 
to  the  circular  cylinder  with  a  diameter /wing  span  ratio  of  0.167. 

The  Reynolds  number  in  the  experiment,  based  on  the  geometric  mean  chord  of  the 
wing,  was  1.0  x  10®.  A  fixed  transition  mechanism  to  stimulate  turbulance  transition 
was  attached  at  12.5  percent  of  the  chord.  The  experimental  data  corresponds  to  a 
free  stream  Mach  number  of  0.4. 

Panelling  is  established  on  only  one  side  of  the  lateral  plane  of  symmetry.  The 
wing  is  modelled  by  11  spanwise  strips,  which  are  located  conveniently  to  match  the 
spanwise  locations  of  the  control  points  with  those  given  in  the  experiment.  Each 
strip  contains  20  upper  and  20  lower  surface  panels.  Axial  panelling  of  the  body 
is  uniformly  spaced  except  between  the  leading  and  trailing  edge  of  the  wing,  where 
panelling  is  matched  to  that  of  the  wing.  The  fuselage  cross  section  panelling  is  equally 
spaced,  with  each  panel  subtanding  an  arc  of  30  degrees.  The  panel  arrangement  of 
the  half  wing-body  configuration  is  shown  in  Figure  4.31. 

The  flow  around  the  wing-body  is  calculated  at  0  and  2  degree  angles  of  attack. 
Neither  compressibility  nor  viscous  corrections  are  made  to  the  calculated  pressures. 

The  calculated  and  experimental  wing  section  pressure  distributions  at  three  dif¬ 
ferent  spanwise  stations  are  presented  in  Figures  4.32  through  4.34.  Effect  of  the 
fixed  transition  is  shown  in  the  pressure  distributions  of  the  experiment  at  12.5  per¬ 
cent  of  the  chord.  The  agreement  is  excellent  in  spite  of  the  fact  that  neither  viscous 
nor  compressibility  corrections  are  applied. 

The  calculated  and  experimental  fuselage  pressure  distributions  are  presented  in 
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Figure  4.30:  RAE  wing-body  configuration. 
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Figure  4.31:  Panel  arrangement  of  the  half  wing-body  configuration. 
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Figure  4.32:  Wing-body  chordwise  pressure  distribution  at  r/R=0.25. 
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Figure  4.33:  Wing-body  chordwise  pressure  distribution  at  r/R=0.60. 
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Figure  4.34;  Wing-body  chordwise  pressure  distribution  at  r/R=0.925. 
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Figures  4.35  through  4.37.  The  present  method  pressure  distribution  is  of  slightly 
lower  magnitude  immediately  above  and  below  the  wing  root.  This  is  believed  to  be 
attributable  to  excursions  from  the  normal  azimuthal  angles  ip  (see  Figure  4.30).  For 
the  present  method,  the  panelling  at  the  wing-fuselage  intersection  is  such  that  the 
actual  values  of  p  are  approximately  ±20  degrees  instead  of  the  nominal  value  of  ±15 
degrees. 

Since  the  panelling  of  the  fuselage  is  matched  to  that  of  the  wing,  which  is  highly 
concentrated  at  the  leading  edge,  axial  lengths  of  the  panels  change  rapidly.  This 
produces  the  steep  variation  of  the  calculated  pressure  near  the  leading  edge.  Again, 
agreement  between  the  calculated  and  experimental  pressure  distributions  is  shown 
to  be  satisfactory. 
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Figure  4.35:  Fuselage  pressure  distribution  in  the  presence  of  wing.  (v?=±15  deg.) 


Figure  4.36:  Fuselage  pressure  distribution  in  the  presence  of  wing.  {<p=±45  deg.) 
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Figure  4.37:  Fuselage  pressure  distribution  in  the  presence  of  wing.  {i^=±75  deg.) 
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Duct  problem 
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4.5  Axisymmetric  duct 


The  next  example  is  chosen  to  assess  the  capability  of  the  present  method  to  predict 
the  internal  flow  properties.  The  geometry  is  a  very  long  axisymmetric  duct  generated 
by  wrapping  a  NACA  0010  section  around  a  circular  cylinder  where  the  cylinder 
length/radius  ratio  is  ten.  The  minimum  internal  cross  sectional  area  is  a  factor  of 
four  smaller  than  the  corresponding  inlet  and  exit  areas. 

Results  for  this  duct  obtained  with  several  panel  codes  were  given  by  Bristow  [2], 
and  additional  results  were  presented  by  Miranda  [19]  and  by  Hess  [lO]  in  a  discussion 
to  Miranda’s  paper.  Predicting  the  pressure  distribution  for  such  an  extreme  duct 
is  a  very  demanding  test  of  a  panel  code.  The  mass  flow  through  the  duct,  which  is 
imposed  by  the  Kutta  condition,  is  extremely  high  in  this  case,  and  there  is  a  tendancy 
for  all  panel  methods  to  underestimate  its  value. 

Given  an  axisymmetric  inflow,  the  solution  is  zixisymmetric,  so  that  the  number  of 
unknowns  is  equal  to  the  niunber  of  chordwise  panels.  The  geometry  is  specified  only 
on  the  meridian  section  in  the  same  way  as  a  two-dimensional  hydrofoil  to  exploit  the 
symmetry,  and  the  influence  functions  are  formed  by  summing  the  individual  panel 
contributions  circumferentially.  Panel  arrangement  for  this  geometry  is  illustrated  in 
Figure  4.38. 

For  this  axisymmetric  case,  the  dipole  strengths  on  the  wake  surface  are  constant 
in  the  circumferential  direction.  The  induced  potential  at  the  control  point  p  on  the 
duct  due  to  the  wake  surface  Sw  can  be  expressed  as 


d  -1 

dxR{p-,q) 


dS 


//  dr  //  11 

(4.1) 

where  Sq  is  a  circular  surface  at  the  exit  of  the  duct  and  S^o  is  a  circular  surface  at 


downstream  infinity.  The  contribution  from  the  surface  at  infinity  is  negligible  as  the 
distance  from  the  control  point  becomes  infinite.  As  a  result,  the  wake  panels  can  be 


113 


Figure  4.38:  Panel  arrangement  for  an  axisymmetric  duct  formed  from  a  NACA  0010 
section  with  a  chord/mean  radius  ratio  of  ten.  (A^c=36,  M#=18) 
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replaced  by  a  circular  dipole  disk  at  the  exit  of  the  duct. 

The  pressure  distribution  for  various  grids  are  shown  in  Figure  4.39.  The  first  four 
calculations  are  for  a  fixed  number  of  chordwise  panels (iVc)  equal  to  36,  but  with  the 
number  of  circumferential  panels  (M#)  varying  between  9  and  60.  The  results  for  36 
and  60  circumferential  panels  are  almost  identical,  and  indicate  a  minimum  pressure 
coefficient  of  -11.3.  Increasing  the  number  of  chordwise  elements  to  60  reduces  the 
minimum  pressure  coefficient  to  -12.5. 

Some  results  obtained  with  other  panel  codes  are  shown  in  Figure  4.40,  which  is 
taken  from  [19].  The  exact  solution,  obtained  with  a  special  high  order  axisymmetric 
panel  code  [3]  shows  a  minimum  pressure  coefficient  of  -13.8  which  is  slightly  lower 
than  the  minimum  value  which  obtained  here.  The  results  obtained  by  a  high  order 
panel  code  developed  by  Hess  [7]  and  by  QUADPAN  [19],  a  low  order  potential  based 
code,  are  similar  to  the  present  method.  The  results  obtained  by  the  Hess  low  order 
velocity  based  code  [6]  show  pressure  minima  much  closer  to  zero,  which  is  evidently 
an  inherent  problem  with  that  method  for  internal  flows. 

An  extension  of  the  present  code  to  the  non-axisymmetric  problem,  with  a  pro¬ 
peller  inside  the  duct  in  a  non-uniform  flow,  is  reported  in  [13]. 
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Figure  4.39:  Pressure  distribution  for  the  axisymmetric  duct  with  the  present  method 
with  different  niimber  of  panels. 
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Figure  4.40:  Pressure  distribution  for  the  axisymmetric  duct  by  different  panel  method 
as  presented  by  Miranda  [19]. 
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4.6  Propeller  performance  analysis  in  steady  flow 

The  final  example  is  a  marine  propeller  on  an  axisymmetric  hub.  The  propeller  is 
NSRDC  4118,  which  is  a  member  of  the  NSRDC  series  propellers  whose  experimental 
results  are  available  from  [l].  The  propeller  is  the  unskewed  propeller  with  three 
blades,  and  the  section  shape  of  the  blade  is  an  NACA  66  mod.  thickness  form 
superposed  on  an  a=0.8  mean  camber  line.  The  geometry  of  the  hub  is  a  circular 
cylinder  with  a  fairwater  whose  ordinates  are  given  as  a  quartic  function  of  the  axial 
distance  from  the  downstream  end  of  the  hub,  as  shown  in  Figure  3.4.  Input  geometric 
parameters  in  Figure  3.4  are  given  as  Rfc=0.2,  X„=0.5,  Xj=0.15,  and  Ar<=0.5. 

Since  the  geometry  and  the  loading  is  repeated  identically  on  each  blade  and  on 
each  inter-blade  segment  of  the  hub,  one  third  of  the  geometry  is  discretized.  Ten 
chordwise  strips  of  panels  are  arranged  to  be  denser  toward  the  tip,  and  each  strip 
contains  twenty  upper  and  twenty  lower  surface  panels. 

Axial  panelling  of  the  hub  is  chosen  to  have  eight  panels  upstream,  twelve  panels 
downstream,  and  twenty  panels  between  the  leading  and  trailing  edges  of  the  blades. 
Eight  circumferential  panels  are  equally  spaced,  with  each  panel  subtending  an  arc  of 
15  degrees.  As  explained  in  Section  3.2,  replacement  of  the  upstream  hub  with  the 
dipole  disk  at  the  upstream  end  of  the  hub  significantly  decreased  the  computational 
axial  distance,  X^,  without  loss  of  accuracy.  Panel  arrangement  of  one  third  of  the 

propeller  is  shown  in  Figure  4.41,  and  that  of  the  whole  propeller  is  shown  in  Figures 
3.2  and  3.3. 

Calculations  are  performed  for  three  different  flow  conditions.  The  design  advance 
coefficient  for  the  propeller  is  J=0.833.  For  the  computations  described  here  the  ap¬ 
propriate  values  of  the  ultimate  wake  radius(R,^),  ultimate  hub  vortex  radius(R,„;.), 
and  tip  vortex  contraction  angle((5c)  were  determined  from  the  experimental  measure¬ 
ments  of  Min  [18].  They  are  given  as  R^=0.83,  Kh=0.1,  and  6,  =  15  degrees.  Distance 
from  the  propeller  plane  to  the  beginning  of  the  ultimate  wake  region(X{,„)  is  given 
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Figure  4.41:  Panel  arrangement  of  one  third  of  the  propeller. 
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as  1.5. 


The  pitch  distribution  of  the  transition  wake  is  given  as  an  input  from  the  results 
by  a  lifting  surface  code  [5],  where  the  wake  is  aligned  to  the  resultant  flow.  Since  the 
ultimate  wake,  which  extends  to  infinity  downstream,  is  replaced  by  the  sink  disk  at 
the  beginning  of  the  ultimate  wake  region,  only  8  panels  between  the  helical  dipole 
strips  are  required  to  represent  the  whole  ultimate  wake.  The  panel  arrangement  for 
the  wake,  including  the  panels  for  the  sink  disk,  is  illustrated  in  Figure  4.42 

The  viscous  drag  coefficient  is  taken  to  be  0.003  for  the  calculation  for  the  design 
condition,  and  that  for  the  off-design  condition  is  taken  from  the  two-dimensional 
section  drag  coefficient  curve. 

The  measured  and  computed  open-water  characteristics  after  the  viscous  correc¬ 
tion  are  shown  in  figure  4.43.  Agreement  is  shown  to  be  satisfactory. 

The  perturbation  potential  distributions  on  the  blade  surface  at  three  different 
radii  at  design  J  are  given  in  Figure  4.44.  The  absissa  is  taken  as  a  nondimensional 
arc-length  along  the  chordwise  strip,  which  is  zero  at  the  lower  trailing  edge  and  one 
at  the  upper  trailing  edge  of  the  blade. 

Blade  pressure  distributions  at  design  J  are  illustrated  in  Figure  4.45  at  three 
different  radii.  The  pressure  distributions  are  shown  to  be  smooth  even  very  close  to 
the  tip.  The  positive  loading  at  the  trailing  edge  near  the  tip,  which  was  experienced 
with  the  other  low  order  panel  method  [15],  is  removed  by  adopting  the  pressure 
Kutta  condition.  The  pressure  difference  between  the  upper  and  lower  surface,  which 
IS  the  loading  on  the  blade,  is  shown  to  have  an  a=0.8  loading  shape  at  this  design 

J.  Unfortunately,  there  are  no  experimental  pressure  data  to  compare  with  these 
computations. 

The  perturbation  potential  distributions  on  the  hub  surface  along  three  different 
streamwise  panels  are  given  in  Figure  4.46.  The  potential  values  at  the  far  upstream 
panels  are  shown  to  be  the  same  in  the  circumferential  direction,  which  is  the  ji.stifi- 
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Figure  4.42:  Panel  arrangement  for  the  wake. 
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Figure  4.43:  Measured  and  calculated  open  water  characteristics  of  NSRDC  propeller 
4118. 
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Figure  4.44:  Perturbation  potential  distribution  on  the  blade  for  NSRDC  propeller 
4118  operating  at  an  advance  coefficient  J=0.833. 
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Figure  4.45:  Computed  chordwise  pressure  distributions  for  NSRDC  propeller  4118 
operating  at  an  advance  coefficient  J=0.833. 
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cation  of  replacing  the  semi-infinite  upstream  hub  by  the  dipole  disk.  The  potential 
jump  at  the  trailing  edge  of  the  root  section  is  preserved  in  the  hub  panels  along  the 
wake  surfzice. 

Computed  hub  pressure  distributions  are  shown  in  Figure  4.47  along  three  different 
streamwise  strips.  Due  to  the  rapid  change  of  the  panel  size  near  the  leading  edge, 
the  pressure  distribution  at  the  leading  edge  is  not  smooth. 
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Figure  4.46:  Perturbation  potential  distribution  on  the  hub  for  NSRDC  propeller 
4118  operating  at  an  advance  coefficient  J=0.833. 
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Figure  4.47:  Computed  hub  pressure  distributions  for  NSRDC  propeller  4118  operat¬ 
ing  at  an  advance  coefficient  J=0.833. 
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Chapter  5 
C  onclusions 


A  surface  panel  method  suitable  for  the  analysis  of  marine  propellers  is  developed  and 
applied  to  various  geometries  to  demonstrate  its  effectiveness.  In  the  way  of  selecting 
the  most  suitable  method  for  the  marine  propllers,  basic  theories  behind  the  various 
panel  methods  are  reviewed  and  characteristics  of  each  method  are  compared.  The 
perturbation  potential  method  is  selected  because  of  its  robustness  to  the  extreme 
geometries,  such  as  a  very  thin  foil  section  and  a  long  duct.  The  perturbation  po¬ 
tential  method  also  requires  relatively  smaller  computing  times  and  computer  storage 
compared  to  the  velocity  method. 

The  usual  Morino’s  Kutta  condition  is  found  to  have  a  deficiency  for  thick  foil 
sections  with  finite  trailing  edge  angles.  Moreover,  it  does  not  account  for  the  three- 
dimensional  cross  flow  effects,  which  become  large  toward  the  blade  tip.  The  deficiency 
for  thick  foil  sections  is  removed  by  including  a  correction  term,  which  is  the  difference 
in  free-stream  potential  values  between  the  trailing  edge  control  points. 

Detailed  study  of  the  flows  at  the  trailing  edge  near  the  tip  suggests  a  pressure 
Kutta  condition,  which  requires  the  pressures  of  the  last  panels  at  the  trailing  edge 
be  equal.  The  non-linear  aspect  of  the  pressure  Kutta  condition  requires  an  iterative 
scheme,  whereby  the  initial  dipole  strengths  on  the  wake  surface  are  obtained  using 
Morino’s  Kutta  condition,  and  the  successive  dipole  strengths  are  adjusted  based 
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on  the  error  in  the  computed  trailing  edge  element  pressures.  Hence,  the  three- 
dimensional  cross  flow  effects  on  Kutta  condition  are  included  by  adjusting  the  dipole 
strengths  on  the  wake  surface  during  the  iteration  process. 

The  compatibility  condition,  which  requires  the  difference  of  the  spanwise  veloci¬ 
ties  of  the  upper  and  lower  panels  at  the  trailing  edge  be  equal  to  the  differentiation  of 
the  spanwise  circulation,  is  shown  to  be  satisfied.  This  indicates  the  self-consistency 
of  the  present  method. 

An  accelerated  iterative  matrix  solver  is  employed  for  the  solution  of  the  resulting 
linear  system  of  equations.  Computing  time  is  reduced  significantly  compared  to  that 
for  the  Gauss  elimination  method. 

Comparison  of  the  results  for  the  ellipsoid  at  zero  angle  of  attack  with  the  analytic 
solutions  provides  a  verification  of  the  present  method. 

Calculation  for  the  circular  wing  with  varying  thickness  is  performed  to  illustrate 
the  thickness  effect  on  the  performance  of  the  wing.  As  the  thickness  is  decreased, 
the  analytic  results  for  the  zero  thickness  wing  by  Jordan  [12]  are  recovered.  This 
demonstrates  the  suitability  of  the  present  method  for  the  analysis  of  a  very  thin 
wing.  With  the  frozen  wake  model  in  the  present  method,  the  thickness  effect  near 
the  tip  is  found  to  have  a  trend  which  is  the  reverse  of  the  two-dimensional  case,  i.e., 
as  the  thickness  is  increased,  the  loading  is  decreased.  This  is  due  to  the  cross  flow 
component  around  the  tip,  whose  magnitude  increases  «is  the  thickness  is  increased. 

The  effect  of  sweep  angles  on  the  spanwise  circulation  of  the  rectangular  wing  is 
investigated.  The  circulation  near  the  tip  is  decreased  for  the  forward  swept  wing, 
while  the  circulation  in  the  inner  part  of  the  wing  is  decreased  for  the  backward  swept 
wing,  compared  to  the  unswept  wing  case.  The  local  section  lift  coefficients  for  the  0 
and  20  degree  swept  wings  show  a  satisfactory  agreement  with  those  by  the  various 
production  codes. 

The  calculated  pressure  distributions  around  the  wing-body  configuration  are  in 
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excellent  agreement  with  those  found  by  experiment  in  spite  of  the  fact  that  neither 
compressibility  nor  viscous  corrections  are  made. 

The  capability  of  the  present  method  to  predict  internal  flow  properties  is  demon¬ 
strated  by  the  example  of  a  long  axisymmetric  duct,  for  which  Hess’s  low  order  surface 
source  method  gives  particularly  inaccurate  results. 

For  the  propeller  problem,  a  rotational  symmetry  is  exploited  so  that  only 
fraction  of  the  geometry,  i.e.,  one  blade  and  the  inter-blade  segment  of  the  hub,  is 
discretized.  By  summing  the  influence  functions  for  the  symmetric  panels  by  rotat¬ 
ing  an  angle  the  same  panelling  on  all  blades  is  achieved.  Substantial  coding 
complication  needed  for  different  panelling  on  the  other  blades  is  thus  avoided. 

A  special  far-wake  approximation  is  employed  for  the  calculation  of  the  influence 
function  of  the  ultimate  wake.  The  semi-infinite  helical  dipole  strip  representing  the 
ultimate  wake  is  replaced  by  a  source  disk  at  the  beginning  of  the  ultimate  wake.  The 
computational  axial  distance  in  the  upstream  hub  is  reduced  without  loss  of  accuracy 
by  replacing  the  semi-infinite  upstream  hub  by  a  dipole  disk. 

A  reliable  pressure  distribution  around  a  marine  propeller,  especially  near  the 
leading  edge  of  the  blade,  is  obtained.  At  the  design  inflow  condition  of  the  selected 
propeller,  the  chordwise  loading  distribution  is  shown  to  have  an  a=0.8  loading  shape. 

The  effect  of  the  hub  is  naturally  included  by  distributing  panels  on  the  surface  of  the 
hub. 

The  viscous  effects  are  included  by  including  a  tangential  friction  force  on  each 
panel  with  the  two-dimensional  section  drag  coefficient.  The  thrust  and  torque,  calcu¬ 
lated  by  the  integration  of  the  element  pressure  on  each  panel,  are  in  close  agreement 
with  the  experimental  results  even  for  the  off-design  conditions. 
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Appendix  A 

Numerical  formulation  of  the 
various  panel  methods  for 
two-dimensional  problems 


To  compare  the  characteristics  of  the  various  methods  described  in  Chapter  2,  each 
method  is  implemented  as  numerical  code  for  the  analysis  of  the  two-dimensional  flow 
around  a  hydrofoil.  Since  the  Hess’s  source  bzised  method  is  widely  used  and  its  char¬ 
acteristics  are  relatively  well  known,  four  other  methods  are  actually  implemented. 


A.l  Perturbation  potential  method 


Two-dimensional  form  of  Equation  2.9,  for  the  field  point  p  on  the  body  surface,  is 


-  TT  (^(p)  =  j 


Sb 


d  d(h 

<Piq)  —  log  R{p]  q)  -  log  R{p;  q) 


ds 


f  d 

+  J  M{q)-^^ogR{p]q)ds. 

Sw  ’ 


(A.l) 


Here  the  three-dimensional  Green  funcion,  is  replaced  by  two-dimensional  one, 

^  log  R,  and  the  surface  integral  is  replaced  by  a  line  integral  along  the  body  and 
wake  surfaces. 
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The  body  geometry  is  replaced  by  an  N-faced  inscribed  polygon,  where  N  is  the 
number  of  panels,  as  shown  in  Figure  A.l.  A  logical  choice  for  the  panel  arrangement 
is  a  cosint  spacing,  where  mean  line  ordinate  and  thickness  is  first  evaluated  at  the 
following  points  along  the  nose-tail  line, 


The  panel  boundaries  are  then  obtained  by  adding  and  subtracting  the  half  thick¬ 


ness  of  the  section  at  right  angles  to  the  mean  line.  This  concentrates  the  elements 
at  the  leading  and  triling  edges,  where  greater  resolution  is  required.  The  nodes  are 
numbered  in  clockwise  order,  starting  at  the  lower  trailing  edge,  and  j-th  panel  is  one 
between  nodes  j  and  j-t-1. 

Singularity  strength  distribution  is  assumed  to  be  piecewise  constant  over  the 
panels.  The  collocation  point,  where  the  discretized  integral  equation  is  satisfied,  is 
selected  as  the  midpoint  of  each  panel.  Then  the  discretized  form  of  Equation  A.l  is 


where 


A  numerical  Kutta  condition  can  be  stated  as  the  potential  jump  in  the  wake 
should  be  equal  to  the  difference  of  total  potential  values,  as  explained  in  Chapter  3. 


{A<f>)yjake  =  =  (i>N  —  4>1  +  Uoo  ^cos  a(x;i^  -  Xi)  +  sin  a{yN  —  r/i)  j , 
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(A.4) 


Figure  A.l 


y 


:  Nomenclature  of  the  potential  methods  for  a  two-dimensional  foil 
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where  C/qq  ^cos  a(i/^  —  ii)  +s\na{yff  —  yi)^  is  the  inflow  potential  difference  between 
the  control  points  of  the  two  trailing  edge  panels. 

This  results  in  a  system  of  linear  equations  for  the  unknown  dipole  strengths. 

-  U^[cosa{xN  -  Xi)  -  sina(y/^  -  Vi)),  i  = 


:=i 


j=i 


(A.5) 


where  =  tt 


,if  *  =  y, 

=  Dij  ,if  *  7^  y, 


and  A{j  =  —  W,-  ,if  i  =  1) 

=  Aij  +  Wi  ,if  I  =  N. 


Solution  of  Equation  A.5  yields  the  values  of  potential  on  the  panels,  under  the 
assumption  that  the  potential  is  constant  on  each  panel.  Surface  velocity  is  obtained 
by  numerical  differentiation  of  the  potential.  A  quadratic  polynomial  to  the  values  of 
potential  at  three  panel  midpoints  is  assumed,  and  the  velocity  at  the  panel  midpoint 
is  obtained  by  differentiating  it  with  respect  to  the  coordinate  that  is  tangent  to  the 
panel.  The  arc-length  between  two  control  points  is  approximated  as  the  sum  of  half 
panel  lengths  of  the  panels. 

Once  the  velocity  is  known,  the  pressure  is  calculated  from  Bernoulli’s  equation, 
where  the  pressure  coefficient  is  defined  as 


P  -  Po 


■ 

Forces  and  moments  are  then  obtained  by  summing  the  element  forces  and  mo¬ 
ments.  An  alternative  lift  coefficient,  which  is  batsed  on  Kutta-Joukowsky’s  law,  is 
calculated  by 


Ct  — 


kfUlc 


=  2 


(A«i) 


wake 


U^c 


and  the  drag  coefficient  should  be  zero. 
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A. 2  Total  potential  method 

Two-dimensional  form  of  Equation  2.12,  for  the  field  point  p  on  the  body  surface,  is 

/Q 

^(9)  —  log  R{p\  <i)ds 

s. 

^^(9)^ — logi2(p;?)ds,  (A.6) 

where  4>oo  =  •  f  =  U^{xcosa  -h  ysino:)  for  uniform  inflow  velocity  with  angle  of 

attack  a. 

Discretizations  of  geometry  and  singularity  strengths  are  identical  to  those  for  the 
perturbation  potential  method.  Then  the  discretized  form  of  Equation  A.6  is 

N 

=  2Uoo{xi  cos  a  4-  y,-  sin  a),  :  =  1,2, . . . ,  N.  (A. 7) 

;  =  l 

A  Kutta  condition  can  be  stated  cls  the  potential  jump  in  the  wake  should  be  equal 
to  the  difference  of  the  total  potential  values. 

(A$),,ait.  =  ^N-  ^1.  (A. 8) 

This  results  in  a  system  of  linear  equations  for  unknown  dipole  strenghs. 

(A.9) 


iV 


^  Aij =  2U^{xi  cos  a  +  y,  sin  a) , 
;=i 


where 

> 

II 

II 

=  A; 

)if  *  7^  J ) 

and 

1 

II 

,if  i  =  1, 

=  Ai,-  +  W, 

,if  =  N 

Here  left  hand  side  of  the  matrix  equation  is  identical  to  that  of  the  perturbation 
potential  method. 
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Calculation  of  the  surface  velocity  follows  the  same  scheme  of  the  perturbation 
method,  except  that  total  velocity  is  calculated  directly  without  summing  the  inflow 
velocity.  The  pressure  and  forces  are  calculated  by  the  same  scheme  of  the  perturba¬ 
tion  method. 

A. 3  Mixed  source  and  vortex  method 


Two-dimensional  form  of  Equation  2.15,  for  the  field  point  p  on  the  body  surface,  is 


TT 


d4> 

dn„ 


^  log -R(p;  9)^5  -  f  — log  R{p]  q)ds - 

J  oug  drip  J  dupdrig 

Sb  Sb 
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Sw 


A4) 


‘■p 


dripdrig 


logJ?(p;  q)ds. 


(A.IO) 


An  equivalent  form  from  Equation  2.16  is 
d<i>  __  f  .  .  d 

Sb 


TT 


dn^ 


d  C  d 

=  j  o{q)— log  R{p;q)ds  +  j  '^{q)  —  [-Q{p-,q)]ds,  (A.ll) 
Sb  ’’  Sb 


where  0(p;  q)  =  arctan(i^)  is  the  two-dimensional  vortex  potential,  <7  =  |^,  7  =  |p, 
and  t  is  defined  to  be  tangential  to  the  body  surface.  Sign  of  the  vorticity  is  defined 
to  be  positive  for  clockwise  induced  flow. 

Equation  A.ll  can  also  be  derived  directly  from  Equation  A.IO  using  the  Cauchy- 

Riemann  relations  for  the  complex  potential  of  log  R  +  iQ, 

d  d 

-logR^—Q, 
ot  on 

Aiogie  =  -|0. 

Since  a  normal  dipole  panel  of  constant  strength  is  equivalent  to  a  pair  of  point 
vortices  at  the  panel  edges,  point  vortices  are  distributed  instead  of  the  piecewise 
constant  dipoles.  Point  sources  are  also  distributed  on  the  same  locations  of  the 


vortices. 
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Locations  of  the  singularities  are  determined  by  first  evaluating  mean  line  ordinates 
and  thickness  at  the  following  points  along  the  nose-tail  line, 

then  adding  and  subtracting  the  half  thickness.  Locations  of  control  points  are  deter¬ 
mined  similarly,  but  x-coordinates  of  the  points  are 

xf  =  -cos(-''^*^  .*  =  l>2,---,Ar/2  +  l. 

In  this  way  the  control  points  are  located  on  the  exact  foil  surface  instead  of  the 
panel  midpoints  which  are  chosen  as  the  control  points  in  the  potential  method. (see 
Figure  A. 2) 

Since  an  explicit  Kutta  condition  is  needed  for  the  velocity  method,  tangential 
velocity  at  the  trailing  edge  control  point  is  set  equal  to  zero.  Direction  of  the  tan¬ 
gential  velocity  is  defined  to  be  normal  to  the  bisector  of  the  upper  and  lower  surfaces 
at  the  trailing  edge. 

Discretization  of  Equation  A. 11  results  in  a  system  of  linear  equations  for  the 
unknown  vortex  strengths. 


where  G,y  =  normal  induced  velocity  at  i-th  control  point 
due  to  the  j-th  point  vortex  of  unit  strength, 
Sij  =  normal  induced  velocity  at  i-th  control  point 
due  to  the  j-th  point  source  of  unit  strength. 


,N, 


(A.12) 


Equation  A.12  is  enough  to  determine  N  unknown  vortex  strengths. 

Surface  perturbation  v'elocity  is  calculated  by  dividing  the  obtained  vortex  strength 
by  the  arclength  between  the  adjacent  control  points.  The  arclength  is  determined  as 

an  arclength  of  the  circle  which  passes  the  vortex  point  and  the  two  adjacent  control 
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y 


Figure  A. 2:  Nomenclature  of  the  velocity  method  for  a  two-dimensional  foil. 
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points.  Surface  velocity  is  obtained  by  adding  the  inflow  velocity  to  the  calculated 
perturbation  velocity.  Calculation  of  the  pressure  follows  the  same  scheme  of  the 
potential  method. 


A. 4  Vortex  method 

Two-dimensional  form  of  Equation  2.18,  for  the  fleld  point  p  on  the  body  surface,  is 


(A.13) 


where  0(p;g)  —  arctan(i^)  is  a  two-dimensional  vortex  potential,  7  =  and  t  is 
defined  to  be  tangential  to  the  body  surface.  Sign  of  the  vorticity  is  defined  to  be 
positive  for  clockwise  induced  flow. 

Discretization  of  the  vortex  method  follows  the  same  scheme  of  the  source  and 
vortex  method.  Instead  of  distributing  both  sources  and  vortices,  only  the  vortices 
are  distributed  on  the  foil  surface. 

Discretization  of  Equation  A.13  results  in  a  system  of  linear  equations  for  the 
unknown  vortex  strengths. 


(A.14) 


where  =  normal  Induced  velocity  at  i-th  control  point  due  to  the  j-th  point 

vortex  of  unit  strength. 

Surface  velocity  at  the  node  point  is  calculated  by  dividing  the  obtained  vortex 
strength  by  the  arclength  between  the  adjacent  control  points.  Calculated  velocity 
here  is  the  total  surface  velocity,  hence  the  inflow  velocity  should  not  be  added.  Cal¬ 
culation  of  pressure  follows  the  same  scheme  of  the  mixed  source  and  vortex  method. 
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Appendix  B 

Equivalence  of  dipole  and  vorticity 
distribution 

Consider  a  surface  S  in  space  bounded  by  a  closed  curve  C,  with  the  unit  normal 
vector  n.  (see  Figure  B.l)  The  vector  between  the  field  point  p(x,y,z)  and  the  source 
point  on  S  is  denoted  R,  and  the  length  of  this  vector  is  denoted  R. 


R=  {x-  +  (y  -  T])j  +  (2  - 

Let  the  surface  S  be  covered  with  a  normal  dipole  distribution  of  strength  then 


there  exists  the  following  relation  between  the  dipole  and  vorticity  distribution. 

THEOREM  :  The  induced  velocity  at  p  due  to  the  normal  dipole  distribution  of 
strength  /z  on  S  is  equal  to  the  sum  of  the  induced  velocities  due  to  a  surface  vorticity 
distribution  on  S  and  due  to  a  line  vortex  along  C.  The  strength  of  the  line  vortex  is 
the  local  value  of  the  dipole  strength,  F  =  ^  (on  C).  The  vorticity  on  S  is  a  vector 
tangent  to  the  curves  of  constant  jj.  and  has  a  magnitude  equal  to  the  surface  gradient 
of  /i  on  S,  -7  =  n,  X  Specifically, 


(B.l) 


c 
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V 


Figure  B.l:  Notation  for  a  general  surface. 
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where  Vp  is  the  gradient  operator  with  respective  to  the  point  p,  and  V,  is  the  gradient 
operator  with  respective  to  the  point  q.  PROOF  :  To  convert  the  line  integral  in 


Equation  B.l  into  a  surface  integral,  start  with  Stokes  theorem, 

jv  -1(11  =  jj{VxV)-ndS. 

c  s 

Put  V  =  ax  P,  where  a  is  an  arbitrary  constant  vector  (a  7^  0),  then 

(a  X  P)  ■  r=  a  •  (P  X  i), 


(B.2) 


and 

V  X  (a  X  P)  •  n  =  n  X  V  ■  (a  X  P)  =  -a  •  (n  X  V)  X  R 
Equation  B.2  becomes 


j  a  - {P  xi)dl  =  -  11  a-{nx'V)x  PdS. 
c  s 

Since  a  is  an  arbitrary  constant  vector,  the  following  transformation  from  a  line  inte¬ 
gral  to  a  surface  integral  is  possible, 

j  t  X  Pdl  =  jj[n  x'V)  X  PdS.  (B.3) 

c  s 

If  P  is  set  to  be  P  =  ^i{q)Vp{^)  in  Equation  B.3,  then 

/fx  dl  =  II {n,  X  V.)  X  dS.  (B.4) 

The  right  hand  side  of  Equation  B.l  is  transformed  into  a  form  which  involves  only 
surface  integrals  on  S,  using  Equation  B.4.  The  integrand  of  the  resulting  surface 
integral  is 


X  n,)  X  Vp(^)  +  (n,  x  V,)  x  Vp(^))  = 


+  V, 
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,  (B.5) 


where  V,  • /iVp(-^)^  means  that  V,  operates  only  on  /iVp(-^),  not  on  n,.  The 
last  two  terms  of  the  right  hand  side  of  Equation  B.5  can  be  expressed  as 

•  '^r(  (5i  ■  '^p(  ^)) . 


and 


n. 


=  ria 


+  n. 


where  V,  •  Vp(-^)  =  ^q('k)  ~  ^  P  ^  Due  to  the  cancellations  of  terms, 
Equation  B.5  becomes 

{B.6) 


li(q)  V.(n,.V,(-i)) 

Equation  B.6  can  be  proceeded  using 

))  =  )  +  n,  X  (V,  X  Vp( 


~R 


(B.7) 


where  the  first  term  is  zero.  Thus  the  integrand  of  the  resulting  surface  integral  in 
the  right  hand  side  of  Equation  B.l  is 


fj.[q)  Hg  X  X  Vp(  ^ 


(B.8) 


On  the  other  hand,  the  integrand  of  the  surface  integral  of  the  left  hand  side  of 
Equation  B.l  is 

f‘(9)V,(n,.V,(^))  =,i{,)n,V,.V,(^)+;<(,)n,x  (V,x  (B.9) 

VVe  can  show  that  Equation  B.9  is  equivalent  to  Equation  B.8,  using 

V,  X  V,(^)  =  -V,  X  (|)  =  -V,i±)  X  i  -  i±)  V,  X  R, 


and 


V,  X  V,(^)  =  V,  X  (|)  =  v,(i)  X  ^  +  (i)  y,  X 


.R31  9V^3 

where  -Vp(^)  =  V,(^),  and  -Vp  xR  =  '7gXR  =  0 
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Appendix  C 

Potential  and  velocities  on  the 
surface  of  an  ellipsoid 


Consider  an  ellipsoid  with  semi-lengths  of  its  axes  a,  b  and  c  in  an  uniform  inflow 
velocity  C/qo  parallel  to  x-axis.  Equation  of  the  ellipsoid  is 


x'  y'  2' 

- h  —  H - =  1. 

a2  62  c2 


(C.l) 


Expression  for  the  perturbation  potential  on  the  ellipsoid  surface  is  given  by  Lamb 


[16], 


where 


d-lx.y.z)  =  U^xK, 


(C.2) 


K  = 
D  = 

Ctn  = 


06c 

2-ao 


D 


d\ 


JO  (a2  +  A)3/3(63  +  A)i/2(c2  +  A)i/2 

abc  D 


The  quantity  D  is  purely  numerical  and  can  be  calculated  by  a  numerical  integration. 
For  brevity,  set  Uoo  =  1,  a  =  6  =  1,  then  Equation  C.l  becomes 

,2 


X  +  y  H - r  —  1. 


(C.3) 
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